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I .  INTRODUCTION 


Time  dependent  two-dimensional  Eulerian  computer  codes  like  HELP1 
2 

and  HULL  are  utilized  to  describe  the  unsteady  interactions  of  continu- 
ous media  (fluid  and/or  solids).  Many  of  these  continuum  codes  have  a 

3 4 5 

common  ancestorial  algorithm,  the  Particle-In-Cell  method  * ’ . During 
the  evolutionary  process,  the  new  codes  have  deviated  substantially  from 
the  original  PIC  method;  for  example,  the  discrete  particles  were  re- 
placed by  a continuum,  certain  Lagrangian-type  features  were  abandoned, 
and  for  calculations  in  solid  mechanics,  material  strength  and  effects 
of  the  deviatoric  stress  tensor  were  included.  These  codes  can  produce 
very  successful  simulations  but  often  the  results  are  not  totally  satis- 
factory. Such  is  the  case  with  the  HELP  code  which  is  used  by  several 
research  laboratories  and  corporations  for  diverse  applications  in  com- 
pressible flow  and  elastic-plastic  flows.  In  certain  ballistics  appli- 
cations at  the  US  Army  Ballistics  Research  Laboratory,  the  code  predicts 
velocities  satisfactorily  but  an  internal  energy  which  implies  a differ- 
ent thermal  state  than  that  indicated  by  experimental  evidence.  The  pur- 
pose of  this  paper  is  to  show  that  the  original  approximations  in  HELP 
lead  to  specific  error  terms  that  significantly  and  consistently  influ- 
ence the  internal  energy  calculation  and  to  propose  a correction  within 
the  context  of  the  current  algorithm. 

The  internal  energy  algorithm  in  HELP  is  based  on  the  finite  differ- 
ence approximation  of  the  total  energy  equation  and  the  kinetic  energy 
calculated  from  updated  mass  and  momentum  values.  This  internal  energy 
approximation  is  shown  to  include  terms  of  the  order  of  the  truncation 
error  which  arise  in  the  kinetic  energy  calculation  from  the  finite  dif- 
ference approximations  of  the  mass  and  momentum  equations.  These  terms, 

T.  Hageman,  L.  J. , et  al . , "HELP,  A Multi-Material  Eulerian  Program 
for  Compressible  Fluid  and  Elastic-Plastic  Flows  in  Two  Space  Dimensions 
and  Time,"  Systems,  Science  and  Software  Report  No.  SSS-R-75-2654, 

July  1975. 

2.  Fry,  M.  A.,  et  al.,  "The  HULL  Hydrodynamics  Computer  Code,"  Air 
Force  Weapons  Laboratory  Report  No.  AFWL-TR-76-193,  September  1976. 

3.  Evans,  M.  W.  and  Harlow,  F.  H.,  "The  Particle-In-Cell  Method  for 
Hydrodynamic  Calculations,"  Los  Alamos  Scientific  Laboratory  Report 
No.  2139,  November  1957. 

4.  Harlow,  F.  H. , "The  Particle-In-Cel 1 Computing  Method  for  Fluid 
Dynamics,"  in  Methods  in  Computational  Physics,  (B.  Alder,  S.  Fernback 
and  M.  Rothenberg,  eds.).  Academic  Press,  New  York,  1964. 

5.  Harlow,  F.  H. , "The  Particle-In-Cell  Method  for  Numerical  Solution 
of  Problems  in  Fluid  Dynamics,"  in  Proceedings  of  Symposia  in  Applied 
Mathematics,  Vol.  XV,  (N.  Metropolis,  J.  Todd,  A.  Tank,  C.  Tompkins, 
eds.)  American  Mathematical  Society,  Providence,  Rhode  Island,  1963. 
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in  certain  cases,  cause  a less  accurate  calculation  of  the  internal 
energy  (they  increase  the  truncation  error)  and  produce  an  interchange 
of  energy  at  each  phase  which  is  not  modeled  in  the  governing  equations. 

Evans  and  Harlow^  identified  an  energy  transfer  mechanism  in  the  convec- 
tion phase  of  the  original  PIC  method  which  can  be  seen  in  HELP.  This 
mechanism  is  due  only  to  spatial  discretization  and  was  illustrated  in 
one  dimension.  The  following  analysis  involves  all  the  phases  in  the 
HELP  algorithm,  is  two-dimensional,  includes  the  effect  of  time  discre- 
tization, and  applies  to  a different  code  with  a different  energy  formu- 
lation (the  PIC  algorithm  transports  total  energy  but  directly  calculates 
the  internal  energy  in  its  other  phases).  Although  this  paper  deals 
exclusively  with  the  HELP  algorithm,  the  concepts  and  results  discussed 
are  applicable  to  other  codes.  In  particular,  the  same  internal  energy 
phenomenon  is  seen  in  calculations  performed  with  the  HULL  code. 

In  Section  II,  the  governing  equations  which  are  modeled  by  the  HELP 
algoritlun  are  listed,  the  corresponding  approximations  are  derived  and 
other  salient  features  of  the  algorithm  are  discussed.  A truncation 
error  analysis  of  the  kinetic  energy  and  the  internal  energy  calculations 
in  Section  III  reveals  the  specific  error  terms  arising  from  the  HELP'S 
finite  difference  approximations  of  the  governing  equations.  In  Section 
IV,  the  modifications  to  the  original  HELP  code  are  given  and  their 
implementation  discussed.  Section  V contains  a mesh  refinement  study 
for  both  the  original  and  modified  version  of  the  code.  A copper  wedge 
impacting  a perfectly  reflective  wall  was  used  for  this  study.  Section 
VI  contains  a detailed  comparison  of  the  two  versions  for  43mm  unconfined 
conical  shaped  charge.  Temperature  profiles  within  the  shaped  charge 
jet  are  given  and  discussed  in  Section  VII.  The  summary  of  the  report 
is  Section  VIII. 


H.  THE  HELP  CODE 

The  unsteady  motion  and  interaction  of  continuous  media  can  be 
described  by  a continuity  equation,  equations  of  motion,  a total  energy 
equation  and  an  equation  of  state.  For  simplicity,  we  shall  consider 
the  Cartesian  formulation.  The  appropriate  two-dimensional  equations  in 
conservative  form  are: 
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where  t,  x,  y,  p,  u,  v,  P,  S , S , S , and  E denote  the  time,  two 

xx  yy  xy 

spatial  coordinates,  density,  x and  y components  of  velocity,  pressure, 
the  two  normal  and  one  shear  stress  components  of  the  stress  deviator 
tensor,  and  specific  total  energy,  respectively.  The  cgs  system  of 
units  is  used  within  the  HELP  code.  However,  in  this  report,  the  re- 
sults are  stated  in  the  SI  system  of  units.  The  elements  of  the  stress 
deviator  tensor  are  functions  of  the  velocity  gradient.  The  pressure 
is  computed  via  an  equation  of  state  of  the  functional  form  P = P(p,I), 
where  I is  the  specific  internal  energy.  The  specific  internal  energy 
is  obtained  as  the  difference  of  the  specific  total  energy  and  the 
specific  kinetic  energy: 


I a E - 0.5  (U2  + V2). 


(S) 


If  the  conservation  eqs.  (1)  - (4)  are  integrated  over  an  arbitrary 
control  area,  the  time  rate  of  change  of  a quantity  within  the  control 
area  can  be  related  to  the  integrals  of  other  quantities  over  the  bound- 
ary enclosing  that  area.  Performing  the  integration  and  using  Green's 
Theorem,  we  obtain 
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where  B is  the  boundary  of  area  A in  the  positive  sense.  Eq.  (7),  for 
example,  equates  the  time  rate  of  change  of  the  x- component  of  the 
momentum  within  the  area  A to  the  product  of  the  specific  momentum  in 
the  x direction  and  the  net  mass  flow  into  the  area  plus  the  sum  of 
certain  surface  forces  (the  pressure  and  the  x-components  of  the  devi- 
ator stress  tensor)  exerted  over  the  boundary  enclosing  the  area  Such 
interpretations  are  used  to  determine  the  HELP  approximations  to  the 
governing  equations. 

The  HELP  code  is  an  Eulerian  code  capable  of  describing  unsteady 
multi-material  interactions  and  of  treating  material  strength  as  an 
elastic-plastic  phenomenon.  A consequence  of  the  multi-material  capa- 
bility is  mixed  cells  (cells  containing  more  than  one  material).  The 
complex  treatment  of  these  cells  is  important  and  indispensable  to  the 
correct  running  of  the  code.  However,  an  accurate  and  complete  analysis 
of  these  numerical  techniques  is  unwieldy.  An  analysis  of  the  pure  cell 
(a  cell  containing  only  one  material)  algorithm  reveals  the  cause  of  the 
internal  energy  problem.  Hence,  the  following  discussion  will  address 
only  the  pure  cell  algorithm.  Furthermore,  we  consider  only  interior 
cells.  We  assume  that  the  grid  spacing  Ax  in  the  x-direction  is  con- 
stant as  well  as  Ay  in  the  y-direction.  The  control  area  A is  taken  to 

be  the  1 , j computational  cell.  See  Fig.  1.  The  left,  right,  top 

and  bottom  boundaries  of  this  cell  are  denoted  by  the  letters  l , r,  a 
and  b,  respectively.  A time  step  At  in  this  explicit  algorithm  is 
determined  by  a Courant  condition.  The  area  integrals  in  eqs.  (6)  - 
(9)  are  approximated  by  m = pAxAy , mu,  mv  and  mF,,  respectively,  where 
m denotes  the  mass  per  unit  length.  All  the  values  are  at  the  center 
of  the  computational  cell.  The  time  derivatives  are  approximated  by  a 
forward  difference.  The  values  of  the  cell  centered  mass,  momentum  and 
specific  total  energy  at  the  new  time  level  are  found  from  the  values 
at  the  previous  time  level.  This  is  accomplished  in  three  stages  by 
determining  the  time  rate  of  change  of  the  mass,  momentum  and  total 
energy  due  to  i)  the  effects  of  the  deviatoric  stresses,  ii)  the 
effects  of  the  pressure,  and  iii)  the  effects  of  the  convection  terms. 
These  phases  are  appropriately  named  SPHASE,  HPHASE  and  TPHASE,  respec- 
tively. During  each  time  step,  each  value  of  the  mass,  momentum  and 
total  energy  is  updated  sequentially  by  each  phase  in  the  order  listed 
and  each  phase  uses  the  previously  updated  values  as  its  initial  values. 
Each  phase  is  solved  independently  of  the  others  and  are  interconnected 
only  through  the  initial  values.  The  boundary  integrals  in  eqs.  (6)  - 
(9)  are  approximated  using  the  current  values  of  the  integrands  at  the 
boundaries  of  the  computational  cell.  The  pressure  is  calculated  before 
SPHASE  and  the  stress  deviator  tensor  is  updated  at  the  beginning  of 
SPHASE.  The  internal  energy  is  updated  by  eq.  (5)  at  the  end  of  each 
phase. 

Specifically,  we  list  the  HELP  approximations  to  the  conservation 
of  mass,  momentum  and  total  energy  equations.  The  HELP  approximations 
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Figure  1.  Computational  Cell  for  Cartesian  Formulation  of  HELP 


to  the  eqs.  (6)  - (9)  in  the  SPHASE  portion  of  the  calculation  for  the 
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where  the  tilde  denotes  the  SPHASE  updated  value  and  the  letter  super- 
scripts r,  H , a,  b refer  to  the  evaluation  of  the  term  at  the  correspond- 
ing boundary.  A variable  without  subscripts  or  superscripts  denotes 

that  quantity  evaluated  at  the  center  of  the  i1^,  jtb  cell  at  the  n1^ 
time  level.  The  boundary  value  of  a variable  is  the  average  of  its  cell 

centered  values  adjacent  to  that  boundary,  for  example,  Sr^x  - 0.3 
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See  l-ig.  1.  The  approximations  (10)  - (13)  can  be  derived  from  a physi- 
cal interpretation  of  the  SPItASr  portions  of  eqs.  (6)  - (9).  For  example, 
consider  the  approximation  (13).  The  effects  of  the  stress  deviator 

tensor  on  the  time  rate  of  change  of  the  total  energy  [mE-mf] /At  during 

the  entire  time  step  At  is  governed  by  the  work  rates  per  unit  surface 

area,  uS  and  vS  acting  on  the  right  and  left  boundar.  ?s  and  vS 

xx  xy  ^ ^ yy 

and  uS  on  the  top  and  the  bottom  boundaries,  times  the  length  of 
these  boundaries. 


The  errors  that  are  introduced  by  the  numerical  approximation  (13) 
can  be  obtained  by  expanding  the  corresponding  finite  difference  equa- 
tion in  a Taylor  series.  Wc  use  the  mass-density  relation  m - p a xa  y 
and  rewrite  cq.  (13)  as 
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A Taylor  series  expansion  of  each  dependent  variable  in  eq.  (14)  at  the 
center  of  the  itb,  jtb  cell  at  the  ntb  time  level  leads  to 
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The  terms  in  eq.  (4)  which  are  relevant  to  SPHASE  are  given  on  the 
left-hand  side  of  eq.  (15)  and  the  dominant  error  terms  appear  on  the 

right-hand  side.  The  order  of  error  terms  are  O(At),  0(Ax2)  and  O(Ay^) . 
Thus,  the  approximation  (13)  is  first  order  in  time  and  second  order 
in  space  within  the  context  of  SPHASE  (assuming  the  post  SPHASE  values 
are  those  at  the  end  of  the  time  step).  A similar  analysis  and  results 
hold  for  approximations  (11)  - (12)  and  eqs.  (2)  - (3),  respectively. 
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The  HPHASE  approximations  are 
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where  the  bar  denotes  the  HPHASE  updated  value  and  P = P(r>,  T).  The 
truncation  error  analysis  of  approximation  (19)  can  be  made  in  exactly 
the  same  manner  as  for  approx  mation  (13).  The  result  is: 
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Thus  eq.  (19)  is  a first  order  approximation  in  time  and  second  order 
in  space  within  the  context  of  HPHASE.  The  other  HPHASE  approximations 
(17)  and  (18)  are  of  the  same  order. 


For  simplicity  in  the  discussion  of  the  TPHASF.  approximations,  we 
assume  that  the  velocity  has  both  positive  x and  y components.  The 
TPHASE  approximations  which  model  the  convection  between  cells  are: 
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where  6m  , 6m  , 5m  , 6m  denote  the  convected  mass  per  unit  length  from 
the  left  and  bottom  cells  and  to  the  right  and  top  cells,  respectively. 

See  Fig.  1.  In  general,  6m^  * p^L^u^At,  where  denotes  the  density 

of  the  cell  from  which  the  mass  is  transported,  is  the  interpolated 

value  of  the  velocity  component  normal  to  the  cell  boundary  and  ld  is 
the  length  of  the  cell  boundary  through  which  the  mass  is  moved.  For 

l l - 1 _ _ 

example,  the  factors  in  6m  are  p » p.  , u = 0.5  (u+u.  , . )/ 

i-l, J i-l,j  . 

[l+At (u  - u.  . ,)/Ax]  and  l/  ■ Ay.  We  note  that  represents  the 

1-  i , J 

transport  velocity  of  6m^  based  on  linear  approximations  over  the  time 

step  At.  The  intuitive  explanation  of  the  TPHASE  approximations  for 

eq.  (21)  is  that  the  mass  at  the  end  of  TPHASE  (the  final  value  at  (n+1) 

time  level)  is  the  mass  originally  in  the  cell  plus  the  mass  transported 

b r 

into  the  cell  (6m  , 6nr)  minus  the  mass  transported  from  the  cell  (6m  , 

6ma) . For  the  total  energy  approximation  ( 24 j , a similar  situation 
exists,  except  that  now  each  convected  mass  is  associated  with  the  speci- 
fic total  energy  of  its  ’’donor"  cell.  Within  the  context  of  TPHASE 

th 

(assuming  the  post  HPHASE  values  are  the  initial  values  at  the  n time 
time  level),  the  approximations  (21)  - (24)  can  be  shown  to  be  first 
order  in  time  and  space  to  the  TPHASE  portions  of  eqs . (1)  - (4),  respec- 
tively. For  example,  to  determine  the  order  of  approximation  (24),  we 
substitute  the  appropriate  mass  approximations  and  obtain 


l 


\ 


A Taylor  series  expansion  of  each  dependent  variable  in  eq.  (25)  at  the 
center  of  the  jr  cell  at  the  nth  time  level  leads  to 


♦ 0(At2)  * 0(Ax2)  ♦ 0(Ay2) . (26) 
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The  terms  in  eq.  (4)  which  are  relevant  to  TPHASE  are  given  on  the 
left-hand  side  of  eq.  (26)  and  the  dominant  error  terms  appear  on  the 
right-hand  side.  The  spatial  gradients  in  the  coefficient  of  At  in  eq. 

-I*  — 0 _ g 

(26)  are  the  results  of  the  transport  velocities'  u , u , v and  v de- 
pendence on  At.  The  order  of  the  error  terms  are  O(At),  O(Ax)  and  O(Ay). 
Thus,  approximation  (24)  is  first  order  in  both  time  and  space,  which 
establishes  the  assertion. 

We  have  shown  that  the  SPHASE  and  HPHASE  approximations  of  eqs.  (1)  - 
(4)  are  first  order  in  time  and  second  order  in  space  and  that  the  TPHASE 
approximations  are  first  order  in  both  time  and  space.  Consequently,  the 
order  of  the  spatial  approximation  in  TPHASE  is  less  than  that  for  either 
SPHASE  or  HPHASE  and  first  order  error  terms  will  be  dominant  within  the 
algorithm. 


III.  THE  KINETIC  ENERGY  AND  INTERNAL  ENERGY  CALCULATIONS 

In  order  to  determine  the  cause  of  the  unphvsical  internal  energy 
values  produced  by  the  1975  HELP  code  in  conical  shaped  charge  calcula- 
tions, we  must  investigate  how  the  total  energy  approximations  are  com- 
bined with  the  kinetic  energy  approximations  to  produce  the  internal 
energy  approximations.  To  this  end  we  list  the  partial  differential 
equations  for  the  kinetic  and  internal  energies.  The  partial  differen- 
tial equation  governing  kinetic  energy  can  be  derived  from  eqs.  (1)  - 
(5)  by  the  following  identity: 
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and  can  be  written  as 
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where  e = 0.5  (u  + v ).  An  interpretation  of  the  above  manipulation  is 
that  given  the  exact  solutions  of  eqs.  (1)  - (5),  the  derived  function  e 
is  identical  to  the  exact  solution  of  eq.  (23).  We  shall  show  that  the 
HELP  approximations  do  not  share  this  property.  The  partial  differen- 
tial equation  governing  internal  energy  can  be  obtained  by  subtracting 
cq.  (28)  from  eq.  (4)  and  by  using  identity  (5): 
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In  the  HELP  code,  the  specific  kinetic  energy  e at  the  end  of  each 
phase  is  computed  via  •••«••* 

e ■ 0,5  j[(mu)/m]2  ♦ [(mv)/m]2^  (30) 

using  the  updated  values  of  the  mass  and  momentum  from  that  phase.  By 
using  the  approximations  (10)  - (12),  (16)  - (18)  and  (21)  - (23),  eq. 
(30)  and  the  expressions  for  the  mass  in  terms  of  the  density,  we  can 
write  the  formulas  used  to  determine  the  updated  specific  kinetic  energy 
at  the  different  phases  in  terms  of  the  initial  values  at  SPHASE,  HPHA3E 
and  TPHASE.  The  resulting  expressions  are  never  explicitly  used  to  cal- 
culate the  kinetic  energy  but  are  numerically  equivalent  to  eq.  (30). 

The  results  are: 
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A truncation  error  analysis  of  eqs.  (31)  - (33)  reveals  significant 
information  about,  the  kinetic  energy  approximations.  Since  we  have 
shown  that  the  dominant  error  terms  within  the  HELP  algorithm  are  of 
first  order,  we  will  not  write  the  higher  order  terms.  Proceeding  in 
a similar  fashion  to  the  truncation  error  analyses  of  Section  II,  we 
obtain  for  eq.  (31) 


+ 0(At“)  ♦ Oi.  Ax“)  + 0 ( Ay~ ) . 


The  terms  in  eq.  (28)  which  are  relevant  to  S PHASE  are  given  on  the  left- 
hand  side  of  eq.  (34)  and  the  dominant  error  terms  appear  on  the  right- 

2 2 

hand  side.  The  order  of  the  error  terms  arc  O(At),  0(/  ”)  and  0(Ay  ). 
Thus,  approximation  (31)  is  first  order  in  time  and  second  order  in 
space,  and  the  order  of  this  approx imat i on  is  consistent  with  the  other 
S PI  IASI:  approximations.  The  terms 
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As  in  SPHASE,  the  HPHASE  kinetic  energy  approximation  is  first  order 
ir.  time  and  second  order  in  space,  and  includes  terms  (those  enclosed  in 
braces  in  eq.  (32))  which  are  of  the  order  of  the  truncation  error.  Eq. 
(36)  shows  that  the  TPHASE  kinetic  energy  approximation  is  first  order 
in  both  time  and  space.  The  terms  enclosed  in  braces  in  eq.  (33)  con- 
tribute only  to  the  0(A:.),  O(Ay),  O(At)  and  higher  order  terms  in  eq. 
(36)  and,  consequently,  are  of  the  order  of  the  truncation  error. 

Thus,  the  order  of  the  approximations  (31)  - (33)  are  in  accord  with 
the  other  approximations  in  the  three  phases  but  these  approximations 
include  terms  which  are  of  the  order  of  the  truncation  error:  terms  of 
order  At  is  SPHASE  and  HPHASE  and  order  At,  Ax  and  Ay  in  TPHASE.  These 
terms  are  consequences  of  calculating  the  kinetic  energy  from  eq.(30) 
and  from  the  particular  choices  made  in  the  finite  difference  approxi- 
mations of  the  mass  and  momentum  equations  in  each  phase.  Furthermore, 
these  terms  do  not  model  any  term  of  the  kinetic  energy  equation.  In 
fact,  if  one  would  write  directly  a finite  difference  approximation  to 
eq.  (28)  in  a consistent  manner  with  the  HELP  approximations  of  eqs.  (1) 
- (4),  the  result  would  be  eqs.  (31)  - (33)  without  the  braced  terms. 
Thus,  the  kinetic  energy  finite  difference  solutions  within  HELP  do  not 
share  the  corresponding  property  possessed  by  the  exact  solution  of  the 
partial  differential  equations:  that  is,  the  function  e,  eq.  (30). 
derived  from  the  finite  difference  solutions  of  the  mass  and  momentum 
equations  does  not  satisfy  the  finite  difference  approximation  of  eq. 
(28).  Although  the  two  approximations  are  the  same  in  the  theoretical 
limit  as  the  mesh  approaches  zero,  in  practice  the  inclusion  of  terms 
of  the  order  of  the  truncation  error  alters  the  accuracy  of  the  calcula- 
tion and  the  computed  value. 

By  casing  the  above  concepts  into  the  framework  of  averaged  quanti- 
ties and  fluctuations  from  their  averages,  an  insigh'  can  be  achieved 
into  the  nature  of  the  truncation  error  terms.  Consider  an  averaging 
procedure  such  that  the  average  of  the  sum  is  the  sum  of  the  averages, 
the  average  of  the  average  is  the  average  and  the  average  of  a fluctua- 
tion is  zero.  The  exact  velocity  can  be  written  as  the  sum  of  the  aver- 
aged velocity  (doubled  barred  quantity)  plus  its  fluctuation  (primed 
quantity).  Componentwise,  we  have 


u = u * u'  and  v * v + v"  . 

The  associated  specific  kinetic  energy  is 

O.S(u2  + v2)  = 0.5(u2  + v2)  * 0. 5(u'2  + u'‘)  ♦ (u  u'  * v u'). 


(37) 
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Using  the  properties  of  the  averaging  procedure,  we  can  write  the  average 
of  eq.  (37)  as 

0.5(u2  + u*)  - 0.5(u2  * v2)  = 0.5(u'~  + v'“).  (38) 
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The  difference  between  the  averaged  exact  specific  kinetic  energy  and 
the  specific  kinetic  energy  of  the  averaged  velocities  is  the  averaged 
specific  kinetic  energy  of  the  fluctuations  which  can  be  called  the 
subgrid-scale  specific  kinetic  energy.  If  we  take  the  averaged  values 
as  the  computed  cell  centered  values  at  the  end  of  a single  time  step, 
then  the  first  term  on  the  left  hand  side  of  eq.  (30)  can  be  associated 
with  the  specific  kinetic  energy  computed  via  the  finite  difference 
approximation  of  eq.  (28)  to  a given  order  of  accuracy  and  the  second 
term  with  the  specific  kinetic  energy  computed  via  the  cell  centered 
values  of  the  mass  and  momentum.  Consider,  for  example,  the  TPHASE 
approximation  to  the  specific  kinetic  energy.  Eq.  (33)  can  be  rewritten 
in  terms  of  the  averaged  values  as 
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^truncation  error  terms| . (39) 
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truncation  error  terms; 


(40) 


Comparing  eq.  (38)  with  eq.  (39),  we  obtain 

{ 

From  eqs.  (40)  and  (39),  we  see  that  the  original  formulation  of  HELP 
excludes  the  subgrid-scale  kinetic  energy.  Accordingly,  the  direct 
calculation  of  the  kinetic  energy  from  eq.  (28)  includes  it. 

The  effect  of  these  truncation  error  terms  is  not  confined  to  the 
kinetic  energy  calculation  but  is  directly  translated  to  the  internal 
energy  calculation  via  eq.  (5).  The  accuracy  of  the  internal  energy 
calculation  is  of  prime  importance,  since  the  pressure,  temperature  and 
strength  properties  of  the  material  directly  depend  on  the  internal  energy 
and  not  on  either  the  total  or  kinetic  energies.  The  dominant  errors  of 
the  internal  energy  approximations  in  SPHASE  can  be  seen  by  subtracting 
eq.  (34)  from  eq.  (IS),  in  HPHASE  by  subtracting  eq.  (35)  from  eq.  (20) 
and  in  TPHASE  by  subtracting  eq.  (36)  from  eq.  (26).  The  results  are 
for  SPHASE, 
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for  TPHASE 


+ 0(At")  * 0(Ax2)  ♦ 0(Ay “)  , 


(42) 


•»  0 (At”)  * 0 (ax”)  * Ofay”) , 


where  wc  have  assumed  Ax  = Ay  in  order  to  simplify  the  TPHASF.  O(At) 
term.  The  accuracy  of  the  internal  energy  calculation  depends  primarily 
on  the  magnitude  of  the  first  order  terms.  The  larger  these  terms  are, 
the  less  accurate  is  the  calculation.  The  structures  of  the  first 
order  terms  in  SPIIASE  and  HPHASE  arc  similar:  a second  time  derivative 
of  (el)  plus  positive  terms  which  are  a consequence  of  the  truncation  er- 
ror terms  in  the  kinetic  energy  calculation.  These  positive  terms  could 
be  excluded  from  the  kinetic  energy  calculation,  and  hence  the  internal 
energy  calculation,  without  changing  the  order  of  the  truncation  error 
of  the  HELP  algorithm.  When  the  second  time  derivative  of  (pi)  is  nor.- 
negative,  the  SPHASF.  and  HPHASK  internal  energy  values  are  computed 
with  less  accuracy  than  would  occur  if  these  truncation  error  terms  were 
excluded.  A typical  time  history  of  the  quantity  (pi)  for  a conical 
shaped  change  simulation  (see  Section  IV)  is  given  in  Fig.  2.  The 
curvature  of  this  graph  is  nonnegative  except  for  a short  interval 
corresponding  to  the  stagnation  region  within  a shaped  charge.  Thus, 
at  least  for  significant  portions  of  a shaped  change  simulation,  the 
internal  energy  calculation  in  SPHASE  and  HPHASE  is  less  accurate 
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figure  2.  Typical  Time  History  Profile  of  the  Internal  1-nergy 
per  Volume  in  a Conical  Shaped  Charge  Calculation 
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because  of  the  inclusion  of  the  truncation  error  terms  from  the  kinetic 
energy  calculation.  Furthermore,  even  if  the  curvature  of  (pI)  with  re- 
spect to  time  was  negative,  the  accuracy  of  the  internal  energy  may  still 
be  less  if  the  order  of  magnitude  of  the  sum  of  the  squared  items  is 
greater  than  t he  curvature  term. 

In  the  TI’llASF  calculation,  each  of  three  first  order  terms  in  eq . 

(.431  should  be  analyzed  to  ascertain  the  effects  of  the  truncation  error 
terms  from  the  kinetic  energy  calculation  (the  squared  terms  in  cq.  (43)) 
on  the  internal  energy  calculation.  In  shaped  charge  simulations,  the  0(/>) 
term  dominates  the  O(Ax)  and  O(At)  terms,  since  the  velocity  and  direc- 
tion of  the  rate  of  change  is  primarily  in  the  y direction  and  since  a 

> 

fairly  coarse  spatial  computing  mesh  must  be  used  ( 0 ( Ay)  = ID  “ versus 
- 8 

0(.At ) -10  ).  Hence,  we  consider  only  the  0(Ay)  term,  which  wc  rewrite 

as : 
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Typical  spatial  profiles  in  shaped  charge  calculations  of  the  quantities 
v and  (el)  along  the  axis  of  symmetry  are  given  in  Fig.  3.  Since  v is 


always  positive  and 
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most  of  their  variation,  the  sum  of  the  first  two  terms  is  generally 
nonnegative.  Consequently,  since  the  third  and  fourth  terms  are  always 
positive,  the  accuracy  of  the  internal  energy  calculation  is  decreased. 
Thus,  as  was  seen  in  the  SPHASF.  and  HPIIASF.  approximations,  t ho  internal 
energy  calculation  for  shaped  charge  calculations  would  generally  be 
more  accurate,  if  the  truncation  error  tcratfrom  the  kinetic  energy  cal- 
culation were  excluded. 


Another  example  of  a flow  in  which  the  truncation  error  terms  in  x Ik- 
kinetic  energy  calculation  would  severely  affect  the  internal  energy 
calculation  is  in  the  expansion  of  perfect  gas  in  a uniform  pressure 

field,  lor  a perfect  gas,  we  have  the  relation  pi  = P(-y-l)  1 where  the 
constant  y is  the  ratio  of  specific  heats.  Consequently,  in  a uniform 

pressure  field,  the  quantities  ? — IilLL.  t lLuL I > UilJJ.  would  be  zero  in 

3t~  ,x 

eqs.  (41)-(43).  If  the  truncation  error  terms  from  the  kinetic  energy 
calculation  were  excluded,  the  entire  0(At)  term  in  SPHASF-  would  be  zero 

as  well  as  the  entire  0(ax)  and  0(Ay)  terms  in  TPIIASF. . Thus  the  internal 

-) 

energy  approximations  in  SPHASF  would  be  0(At")  and  TPHASF.  0(ax“)  and 

0 (A y “ ) . Consequently,  the  approximations  achieve  a higher  order  of 
accuracy  when  the  extraneous  terms  are  excluded. 
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Figure  3.  Typical  Spatial  Profiles  of  the  Axial  Velocity  anti  Internal 
Energy  per  Volume  in  a Conical  Shaped  Charge  Calculation 
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In  other  applications,  the  order  of  magnitude  and/or  the  algebraic 
signs  of  the  first  order  terms  in  eqs.  (41)  - (43)  must  be  analyzed  in 
order  to  determine  the  effects  of  truncation  error  terms  in  the  kinetic 
energy  calculation  on  the  accuracy  of  the  internal  energy  calculation. 

In  regions  of  large  gradients,  the  magnitude  of  these  truncation  error 
terms  can  be  large  because  of  their  quadratic  dependence  on  the  first 
spatial  derivatives.  Since  limitations  on  running  time  and  machine  stor- 
age necessitate  fairly  coarse  computing  meshes  for  two-dimensional  simu- 
lations, the  truncation  errors  related  to  the  finite  size  of  the  mesh 
cell  are  more  likely  to  be  important  than  those  related  to  the  time  step. 
Thus,  the  0(ax)  and  O(Ay)  terms  in  eq.  (43)  may  dominate  the  truncation 
error.  We  shall  see  that  this  is  the  case  in  shaped  charge  calculations. 

The  inclusion  of  the  truncation  error  terms  in  the  kinetic  energy 
calculation  alters  not  only  the  computed  values  of  the  kinetic  energy  of 
a cell  at  each  cycle  but  also  the  values  of  the  internal  energy.  The 
coefficients  of  the  At  terms  in  eqs.  (31)  and  (32)  are  positive  and  in- 
crease the  kinetic  energy.  For  equal  spatial  meshes  (Ax  = Ay),  the  entire 
first  order  term  in  the  TPHASE  calculation  is  negative  for  Courant  num- 
bers less  than  a half  and  decreases  the  kinetic  energy.  The  effect  of 
truncation  error  terms  on  the  internal  energy  is  reversed  because  of  eq . 
(5) . The  SPHA5E  and  HPHASE  terms  decrease  the  internal  energy  and  the 
TPIIASE  increases  it.  Thus,  these  terms  can  be  interpreted  as  a transfer 
mechanism  which  is  not  modeled  by  the  governing  equations  and  which  con- 
verts internal  energy  into  kinetic  energy  and  kinetic  energy  into  internal 
energy.  Consider,  for  example,  the  one  dimensional  first  order  energy 
approximation  in  TPHASF.  for  motion  in  the  x-direction.  The  only  first 
order  term  that  the  interna]  energy  calculation  includes  is  the  posi- 
tive term: 


•Vi  u 

:pnM 


(A 


( 45) 


Expanding  expression  (45)  in  a Taylor  series  about  the  cell  center  and 
the  n**1  time  level,  we  obtain 


(461 


i 2 

j . to  the  lowest  order,  where  X = 0.5puAx  and  X'  = <-0.5pu  At.  If  (X  - X') 

: were  the  coefficient  of  viscosity,  then  expression  (46)  would  be  identi- 

t cal  to  the  viscosity  term  in  the  one-dimensiona’  internal  energy  cqua- 

[ tion  for  a viscous  fluid.  Thus,  the  energy  transfer  mechanism  in  this 

case  could  be  defined  as  an  explicit  artificial  viscosity  term,  since  it 
is  explicitly  included  in  the  difference  equation  much  like  the  implemen- 
tation of  the  von  Neumann  and  Richtmyer  artificial  dissipation  scheme. 
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Evans  and  Harlow  identified  the  term  corresponding  to  >.(3u/3x)  in 

their  one-dimensional  analysis  of  the  original  PIC  code.  The  v cm 

2 

* (3u/3x)  is  not  included  in  their  analysis,  since  they  did  not  include 

the  effect  of  time  differencing.  From  expression  (46)  , we  see  that  the 
time  discretization  decreases  the  amount  of  kinetic  energy  converted  to 
internal  energy  in  TPHASE.  This  explicit  type  of  artificial  viscosity 
is  confined  only  to  the  TPHASE  energy  calculation  and  is,  in  addition 

to  the  implicit  artificial  viscosity6  (that  type  of  artificial  viscosity 
deduced  from  the  neglected  truncation  error  terms  within  an  algorithm), 
already  inherent  in  a first  order  algorithm. 


I V . THE  MODIFIED  HELP  CODE 

We  have  shown  that  the  terms  of  the  order  of  the  truncation 
error  which  are  included  in  the  kinetic  energy  calculation  decrease 
the  accuracy  of  the  internal  energy  calculation  under  certain  circum- 
stances. To  determine  the  effects  of  omitting  these  terms  in  such  a 
calculation,  the  HELP  code  was  modified  to  allow  a kinetic  energy  calcu- 
lation which  did  not  include  the  first  order  terms  in  eqs.  (31)  - (33). 
Consequently,  the  kinetic  energy  was  not  computed  by  the  updated  mass 
and  momentum  values  but  was  considered  a separate  dependent  variable. 
This  kinetic  energy  for  both  pure  and  mixed  cells  was  updated  according 
to  a direct  finite  differencing  of  eq.  (28)  in  a manner  consistent  with 
the  other  approximations  and  was  stored  in  an  array.  The  array  TKF.G 
contains  the  cell  centered  specific  kinetic  energy  and  the  array  TKF.GM 
contains  the  specific  kinetic  energy  of  each  material  in  an  interface 
cell.  By  implementing  such  a scheme,  the  algorithm  still  numerically 
conserves  mass,  momentum  and  specific  total  energy  as  before  but  now 
possesses  a new  discrete  dependent  variable, the  specific  kinetic  energy, 
and  another  finite  difference  equation.  However,  the  modified  version 
required  slightly  less  computing  time  than  the  original  version,  since 
all  the  quantities  needed  to  compute  a new  kinetic  energy  value  are 
already  available  from  the  mass  and  momentum  calculations  and  future 
references  to  the  kinetic  energy  are  simply  retrievals.  Throughout 
the  modified  algorithm  refert.nces  to  the  specific  kinetic  energy 
are  always  to  the  arrays  TKEG  and  TKEGM  and  never  to  the  quantity 
2 2 

0.5  (u  + v ).  Although  other  formulations  are  possible,  the  present 
one  calculates  the  specific  kinetic  energy  in  the  desired  manner  (hence, 
the  internal  energy)  and  is  relatively  easy  to  implement  in  the  existing 
HELP  code. 

Since  an  accurate  calculation  of  the  specific  internal  energy  is  the 
ultimate  goal,  a natural  alternative  to  the  1975  HELP  algorithm  is  the 


6.  Roachc,  P.  .1.,  Computational  Fluid  Dynamics,  Chapter  V,  Hermosa, 
Albuquerque,  1976. 
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direct  calculation  of  the  specific  internal  energy  in  every  phase.  This 
type  of  algorithm  deletes  the  truncation  error  terms  ir.  the  internal  er- 
ror calculation.  In  several  examples  (shaped  charge  simulations  were  not 
7 

included),  Karpp  obtained  significantly  better  internal  energy  values 
with  this  approach.  A closely  related  algorithm  is  the  1971  version  of 

g 

HELP.  This  version  directly  solves  for  the  specific  internal  energy 
in  the  SPHASE  and  HPHASE  portions  of  the  algorithm  but  transported  the 
specific  total  energy  in  TPHASE.  A stated  reason  (see  Ref.  1)  for  the 
current  total  energy  version  of  HELP  was  the  problems  rising  from  the 
internal  energy  calculation  in  HPHASE.  The  two  pass  energy  calculation 
in  HPHASE  caused  inaccurate  definitions  of  pressure  and  velocities  near 
free  surfaces  and,  consequently,  the  inclusion  of  the  GLUE  subroutine 
to  rectify  them.  The  GLUE  subroutine  was  not  totally  satisfactory.  The 
desire  to  avoid  these  past  difficulties  compelled  us  away  from  an 
internal  energy  formulation  and  to  the  retention  of  the  total  energy 
formulation. 

The  modifications  listed  below  apply  to  the  version  of  HELP  described 
in  Ref.  1,  This  original  formulation  is  available  on  the  RRL  GDC  7600 
under  the  permanent  file  name  BRLHELP,  CY  = 14.  The  permanent  file 
identification  is  ID  = CMCHYDRG.  The  modification  deck  is  listed  in 

9 

Appendix  A.  It  is  in  CDC  UPDATE  format  and  contains  comment  cards  in 
order  to  explain  each  particular  change.  The  modified  version  of  HELP 
is  available  on  the  BRL  CDC  7^00  under  the  permanent  file  name  BRL11EU', 

CY  = 15.  The  permanent  file  identification  is  ID  = CMCUYDRO. 

In  addition  to  the  changes  in  the  common  deck,  IIFIPCOM,  the  modi- 
fications affect  15  of  the  57  subprograms  in  HELP:  namely,  SIM  IASI. , 
HPHASE,  TPHASE,  EDIT,  EQST,  INPUT,  SETUP,  SETUPA,  Fll.GRD, FLGSET, 

RNDOFF , NEWMIX,  NEWFLG,  ENCHCK,  and  UVMOD.  The  modifications  incorpo- 
rated into  HELP  enables  one  to  run  simulations  involving  all  the  options 
of  HELP  described  in  Ref.  1 except  the  penetration  and  rezonc  packages. 
The  modifications  have  not  been  included  in  any  of  the  plugging  rou- 
tines: PLGADD,  PLCALF,  PLGGEN,  PLGMAS,  PLGTCR,  PLGVOL,  Pl.UGtlV,  PTSAV, 


7.  Karpp,  R.  R.  and  Soldstein,  S.,  "Modifications  to  the  Hydrodynamic 
Computer  Code  HELP,"  Los  Alamos  Scientific  Laboratory  Summary  Report 
R199,  M-4-1S96,  December  1977. 

8.  Hagcman,  L.  J.  and  Walsh,  .1.  M. , "HELP,  A Mult i -Material  Eulerian 
Program  for  Compressible  Fluid  and  Eiast ic-Plast ic  Flows  in  Two  Space 
Dimensions  and  Time,  Vol . I,  Ballistic  Research  Laboratory  Contract 
Report  NO.  39,  May  1971.  (AD  #726459) 

9.  Control  Data  Corporation,  "UPDATE  Reference  Manual,"  Control  Data 
Corporation  Manual  60-449900,  1975. 


or  in  the  REZONE  subroutine.  However,  a rezoning  deck 
in  conjunction  with  the  REZONE  subroutine  was  obtained 
and  subroutine  were  modified  to  allow  a limited  rezone 
the  modified  formulation.  The  restrictions  on  this  re 
i)  the  rezone  package  must  be  user-activated,  ii)  the 
grid  can  be  extended  only  in  the  axial  direction,  and 
row  of  computation  cells  (J  = JMAX)  must  contain  void, 
package  combines  cells  as  described  in  the  HELP  manual 
Jock  tIDEN'Ts  REZJAS  and  CHREZ)  and  the  modifications  ( 
are  listed  in  Appendix  B in  CDC  UPDATE  format. 


which  is  used 
. Both  the  deck 
capability  in 
zone  package  are: 
computational 
iii)  the  last 
The  rezone 
The  re zone 
INDENT  CORREZ) 


The  modifications  are  implemented  in  such  a way  as  to  retain  the 
total  energy  conservation  within  the  computational  grid.  At  the  end 
of  859  cycles  (t  = 16  gs)  in  the  43  mm  shaped  charge  calculation,  pre- 
sented in  Section  VI,  the  maximum  relative  error  between  the  theoretical 

.9 

total  energy  and  the  actual  total  energy  within  the  grid  is  1.09  x 10 
lor  this  calculation,  the  modifications  increased  the  storage  require- 
ments by  only  6*0  and  decreased  the  running  time  by  approximately  2%. 

Finally,  we  note  that  another  correction  set  is  listed  in  Appendix 
A under  t he  TDF.NT  CORMAP.  This  enables  the  MAPS  subrout:-.e  to  indicate 
small  magnitude  negative  values  in  the  MAPS  output.  Pr  r to  this 
correction,  a single  dash  was  used  as  a symbol  for  both  positive  and 
negative  ranges  of  a variable's  value.  With  this  correction  a single 
dash  represents  positive  values  and  a double  dash  negative  values. 


V.  MESH  REFINEMENT  STUDY  FOR  A WEDGE  IMPACT  CALCULATION 

Since  the  study  of  the  internal  energy  problem  is  based  on  a trunca- 
tion error  analysis,  a mesh  refinement  study  is  appr  -iate  to  the  inves- 
tigation of  the  effects  of  the  truncation  error  terms  on  the  accuracy  of 
the  computed  specific  internal  energies.  Because  of  the  complexity  of 

a full  shaped  charge  simulation,  a related  model  problem1  is  simulated 
in  the  mesh  refinement  study.  An  observer  stationed  at  the  stagnation 
point  in  Fig.  7b  (page  37)  would  find  the  uncollapsed  liner  moving 
towards  him  and  separating  into  two  parts  (the  slug  and  jet).  Conse- 
quently, a 20  mm  wide  copper  wedge  traveling  at  2 km/s  was  simulated 
as  it  impinged  on  a perfectly  reflective  wall  at  30°  of  obliquity.  The 
wedge  collapse  was  modeled  with  slab  symmetry.  The  material  constants 
for  the  copper  arc  CZERO  = 0.235  C.Pa,  RMl)  = 45.50  CPa,  STK1  = 0.95  CPa, 


10.  l.acetera,  .1.  M. , USA  Ballistic  Research  Laboratory,  Private 
Communication,  April  1978. 

11.  Rirkhoff,  C.  , et  al.,  "Explosives  with  Eined  Cavities,"  Journal 
of  App 1 i cd  Physics  19,  563-582,  1948. 
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STK2  =5.5  GPa,  STE2  = 0.53  MJ/kg  and  AMDM  = 0.9785.  See  the  HELP 

manual*  for  the  definition  of  the  above  constants.  The  initial  density 

of  the  copper  is  8.9  x 10^  kg/m^.  We  note  that  the  initial  pressure  and 
specific  internal  energies  are  zero,  since  the  code  computes  the  over- 
pressures and  specific  internal  energies  above  the  initial  state.  Con- 
sequently, throughout  this  report,  the  values  of  the  specific  internal 
energy  are  relative  to  the  initial  reference  value  and  are  not  the 
absolute  values.  The  material  surrounding  the  wedge  is  void.  A 

1 2 

single  pass  through  SPHASE  was  used.  The  Tillotson  equation  of  state  c 
for  copper  was  used.  The  coarse  computational  mesh  was  50  by  100  cells 
(Ax  = Ay  = 1 mm)  and  the  fine  mesh  was  60  by  200  cells  (ax  = Ay  = 0.5 
mm).  Both  the  original  and  modified  version  of  the  HELP  code  were  run 
on  r>  C,DC  7600. 

The  results  at  10us  after  the  wedge  first  impacts  the  wall  arc  shown 
in  Figs.  4-6.  Besides  the  specific  internal  energy,  two  other  important 
quantities,  the  relative  axial  velocity  (relative  to  the  end  of  the  slug 
portion)  and  density,  are  compared  along  the  wall  in  the  slug,  stagna- 
tion and  jet  regions.  Qualitatively,  the  entire  compression  and  rela- 
tive axial  velocity  curves  and  the  slug  port;on  of  the  specific  internal 
energy  curve  are  similar  among  the  mesh  variations  and  version  changes. 
However,  the  qualitative  behavior  of  the  jet's  specific  internal  energy 
markedly  differs  between  the  two  versions  for  either  mesh.  This  is  due 
to  the  increased  gradients  within  the  jet  portion  which  can  drastically 
increase  the  magnitude  of  the  truncation  error  terms  from  the  kinetic 
energy  calculation,  and  thus,  the  truncation  error.  Quantitative  com 
parisons  were  made  at  three  positions  in  Figs.  4-6:  the  slug  end,  the 
stagnation  point  and  the  jet  end.  The  symbol  4 denotes  the  stagnation 
point  (the  cell  along  the  wall  with  the  largest  pressure.  While  the 
total  deviation  among  the  relative  axial  velocity  and  compression  curves 
from  the  largest  computed  value  at  the  three  stations  was  less  than  71 
(except  for  9.91  in  the  compression  at  the  jet's  end),  the  total  devia- 
tion of  the  specific  internal  energy  curve  was  82.11,  26.41  and  651  at 
the  jet,  stagnation  and  slug,  respectively.  The  tremendous  increase  in 
the  total  deviation  of  the  specific  internal  energy  curve  is  due  to  the 
first  order  terms  included  in  the  calculation  of  the  kinetic  energy.  In 
fact,  the  deviation  of  the  coarse  mesh  specific  internal  results  from  the 
fine  mesh  results,  decreased  by  over  t>2l  from  the  original  version  to  the 
modified  version.  Consequently,  one  obtains  much  less  variation  in  the 
specific  internal  energy  values  as  one  converges  to  the  solution  through 
a mesh  refinement  with  the  modified  version.  Fig.  6 suggests  that  the 
specific  internal  energy  values  which  are  computed  by  the  modified  ver- 
sion and  with  the  coarse  mesh  provide  a better  approximation  to  the  exact 
values  than  those  values  computed  by  the  original  version  with  the  fine 


12.  Tillotson,  J.  II.,  "Metallic  Equation  of  State  for  llypcrvo  loe  i t y 
Impact,"  General  Atomic  Report  No.  GA-5216,  July  1962. 


RELATIVE  AXIAL  VELOCITY  (km/s) 


Figure  4.  Spatial  Profiles  of  the  Relative  Axial  Velocity  at  lOys 
Along  the  Wall  for  a Copper  Wedge  Impact  Calculation 
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Figure  5.  Spatial  Profiles  of  the  Compression  at  10gs  Along  the 
Wall  for  a Copper  Wedge  Impact  Calculation 


mesh.  This  trend  is  also  present  in  the  compression  curve.  The  results 
of  Fig.  6 also  show  that  for  this  type  of  problem,  the  truncation  error 
terms  in  the  kinetic  energy  calculation  associated  with  TPHASE  (Courant 
number  of  0.4  is  used)  dominate  those  associated  with  SPHASE  and  HPHASE 
and  cause  the  specific  internal  energy  to  be  increased. 


VI.  THE  COMPARISON  OF  THE  MODIFIED  AND  ORIGINAL  FORMULATIONS 
“"USING  A 43MM  CONICAL  SHAPED  CHARGE  SIMULATION 


We  will  consider  the  43mm  conical  shaped-charge  which  is  unconfined  and 
has  a rounded  apex.  See  Fig.  7a.  The  actual  warhead  is  obtained  by 
rotating  Fig.  7a  about  the  axis  of  symmetry.  The  explosive  is  detonated 
and  the  detonation  wave  collapses  the  conical  liner  towards  the  axis  of 
symmetry  with  a varying  velocity.  Sixteen  microseconds  after  detonation, 
the  liner  consists  of  three  parts  (Fig.  7b'':  The  uncollapsed  liner,  the 
low  velocity  large  mass  slug  and  the  high  elocity  small  mass  jet.  By 
performing  a Galilean  transformation  at  the  point  where  an  original  ring 
of  liner  impinges  on  the  axis  of  symmetry,  a stagnation  point  in  the 
flow  can  be  shown  to  exist.  This  stagnation  point  divides  the  collapsed 
liner  into  the  slug  and  jet.  It  is  the  jet  which  pierces  the  armor  and 
is  of  prime  concern  to  the  shaped  charge  designer.  By  assuming  that  the 
liner  can  be  modeled  as  an  incompressible  material  under  typical  loading 

conditions,  several  researchers^ ’ have  developed  one-dimensional 
analytical  or  simple  numerical  models  to  determine  the  velocity  field. 
However,  when  two-dimensional  axisymmetric  geometry,  compressibility, 
strength,  and  thermodynamic  effects  are  included,  a full  numerical  simu- 
lation is  necessary.  The  liner  collapse  and  early  jet  formation  of  this 

shaped  charge  was  simulated  using  the  BRLSC15  code  (the  precursor  of  the 
1975  HELP  code)  and  the  results  compared  to  experiments  in  Ref.  16. 


13.  Pugh,  E.  M. , Eichelberger,  R.  M. , and  Rostoker,  N. , "Theory  of  Jet 
Formation  by  Charge  with  Lined  Conical  Cavities,"  Journal  of  Applied 
Physics  19,  563-582,  1948. 

14.  Kiwan,  A.  R.  and  Wisniewski,  H.,  "Theory  and  Computations  of 
Collapse  and  Jet  Velocities  of  Metallic  Shaped  Charge  Liners,"  US  Army 
Ballistic  Research  Laboratory  Report  No.  1620,  1972.  (AP  #907161L) 

15.  Gitting,  M.  L.  , "BRLSC:  An  Advanced  F.ulerian  Code  for  Predicting 
Shaped  Charges,  Volume  1,"  US  Army  Ballistic  Research  Laboratories 
Contract  Report  No.  279,  December  1975.  (AD  #A023962) 

16.  Harrison,  J.  T.,  "A  Comparison  Between  the  Equlerian  Hydrodynamic 
Computer  Code  (BRLSC)  and  Experimental  Collapse  for  a Shaped  Charge 
Liner,"  US  Army  Ballistic  Research  Laboratory  Memorandum  Report  No. 
ARBAL-MR- 02841,  June  1978.  (AD  #A059711) 
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Figure  7a.  Initial  Conical  Shaped  Charge  Configuration 
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Figure  7b.  Copper  Configuration  at  16ns  for  Conical  Shaped  Charge 
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Since  the  HELP  code  is  now  generally  used  instead  of  the  BRLSC  code  for 
shaped  charge  studies,  we  shall  limit  our  discussion  to  the  HELP’S  re- 
sults. In  order  to  further  compare  the  original  and  modified  versions 
of  the  HELP  code,  a conical  shaped  charge  with  a copper  liner  and  Compo- 
sition B explosive  (Pig.  7a)  is  modeled  in  cylindrical  coordinates  and 
with  identical  material  constants,  initial  conditions,  code  options  and 
grid  structure.  The  material  constants  for  the  copper  are  identical  to 
those  used  in  the  copper  wedge  calculation.  The  slipline  option  is 
employed  along  the  copper  - Composition  B interface.  A void  surrounds 
the  entire  warhead.  The  initial  state  is  quiescent:  all  properties 

3 3 

zero  except  the  standard  values  of  the  densities  (8.9  x 10  kg/m')  for 

3 3 

copper  and  1.717  x 10  kg/m  for  Composition  B.  The  explosive  is  point 
detonated  at  the  intersection  of  the  axis  of  symmetry  and  the  base  of 
the  Composition  B.  The  computational  mesh  is  60  cells  in  the  radial 
direction  and  187  cells  in  the  axial  direction.  The  axial  increment 
is  A z = 0.52  mm  and  the  radial  increment  is  Ar  = 0.52  mm  for  the  first 
46  cells  and  slowly  increases  to  19.93  mm  at  the  last  radial  cell.  In 
both  calculations  the  bottom  boundary  is  transmitt ive.  The  strength 
option  with  a single  pass  is  utilized  (that  is,  the  contribution  of  the 
stress  deviator  tensor  is  included).  The  Jones-Wi lkins-Lee  equation 

17  12 

of  state  for  the  explosive  and  the  Tillotson  equation  of  state  for 

the  copper  are  used.  The  constants  used  in  these  equations  of  state 

are  given  in  the  HELP  reference  manaual*.  Sample  input -output  for  this 
shaped  charge  is  given  in  Ref.  18.  The  calculations  were  performed  on 
a CDC  7600. 

Fig.  8 shows  the  specific  kinetic  energy  along  the  axis  of  symmetry 
from  the  slug  tail  to  the  jet  tip  for  both  formulations  at  1 6v s . The 
symbol  + denotes  the  stagnation  point.  This  comparison  shows  that  modi- 
fied algorithm  values  are  always  larger  than  those  of  the  original 
algorithm.  Thus,  the  truncation  error  terms  in  the  kinetic  energy  calcu- 
lation associated  with  the  spatial  discretization  dominate  those  associated 
with  the  temporal  discretization.  In  fact,  throughout  the  flow  field, 
the  original  formulation  gave  smaller  values  of  kinetic  energy  than  the 
modified.  Fig.  9 is  the  comparison  of  the  specific  internal  energies 
along  the  axis  of  symmetry  at  16us.  The  modified  code  predicts  up  to 
88%  decrease  in  the  specific  internal  energy  values  of  the  original 
code.  Along  the  axis,  the  original  formulation  predicts  a liquid-vapor 
slug  and  a liquid-vapor  jet,  whereas,  the  modified  formulation  predicts 


17.  Lee , E . , Finger,  M.  and  Collins,  W. , "JWL  Equation  of  State 
Coefficients  for  High  Explosives,"  Lawrence  Livermore  Laboratory 
Report  No.  UCID-16189,  January  1973. 

18.  Lacetera,  J. , Jr.,  Lacetera,  J.  M.  and  Schmitt,  J.  A.,  "The  BRL 
7600  Version  of  the  HELP  Code,"  USA  Ballistic  Research  Laboratory 
Report  (in  preparation). 
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of  the  Specific  Internal  Energies  Along  the  Axis  of  Symmetry  at  16us 
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a solid  slug  and  a jet  which  remains  a solid  until  near  the  tip  region 

where  it  liquifies.  We  note  that  the  BRLSC  code  results16  indicate 
even  a higher  specific  internal  energy  value  than  the  HELP  1975  algo- 
rithm. rhe  specific  internal  energy  value  corresponding  to  the  melting 
point  of  coppe?  53  MJ/kg)  is  a material  property  input  value.  The 
specific  intern  energy  value  corresponding  to  the  incipient  vaporiza- 
tion point  of  copper  (1.38  MJ/kg)  is  contained  within  the  Tillotson 
equation  of  state.  In  the  original  formulation,  the  specific  internal 
energy  of  each  cell  in  the  entire  jet  is  well  above  the  value  corres- 
ponding to  incipient  vaporization  at  16us.  Thus,  the  jet  is  character- 
ized as  a liquid-vapor  jet.  However,  in  the  modified  formulation,  the 
specific  internal  energy  of  the  same  cells  is  generally  below  that  for 
melt  and  a solid  jet  with  several  melted  sections  is  predicted.  See  Fig. 
10.  Although  no  actual  temperature  measurements  have  been  taken  for 

19 

this  shaped  charge,  experimental  evidence  of  Von  Holle  and  Trimble 
supports  the  conclusion  that  the  copper  jet  is  in  the  solid  state. 

Thus,  qualitative  thermal  agreement  is  achieved  for  the  first  time  with 
the  modified  formulation.  Fig.  11  illustrates  the  specific  internal 
energy  values  generated  by  the  modified  formulation  at  various  times 
along  the  axis  of  symmetry  from  just  behind  the  stagnation  point  to  the 
jet  tip.  The  relative  positions  of  the  curves  do  not  indicate  their 
actual  relative  positions  predicted  by  the  code.  The  decrease  in  the 
specific  internal  energy  values  in  the  jet  as  a function  of  time  is  not 
constant  and  become  less  as  time  increases.  The  proportion  of  the  com- 
putational cells  within  the  jet  with  a specific  internal  energy  value 
higher  than  that  of  melt  decreases  from  1.00  to  0.76  to  0.32  to  0.24  as 
time  increases  from  8 to  12  to  16  to  18.5us,  respectively. 

Major  differences  in  the  internal  energy  illustrated  in  Fig.  9 will 
drastically  affect  the  material  properties  of  the  slug  and  jet.  An 
investigation  into  the  copper's  strength,  ductibility  and  cohesion 
cannot  be  computationally  investigated  with  the  original  version  of 
HELP,  since  it  predicts  a liquid-vapor  state  for  most  of  the  copper  in 
the  slug  and  jet.  In  particular,  when  the  specific  internal  energy  of 
a cell  is  equal  to  or  greater  than  chat  for  melt,  the  stress  components 
are  equated  to  zero  and  the  flow  variables  remain  unchanged  during 
SPHASE.  Consequently,  in  the  original  formulation  and  for  the  43mm 
conical  shaped-charge  problem,  the  inclusion  of  the  strength  option 
does  not  generate  significantly  different  values  of  the  flow  variables. 
Since  the  values  of  the  specific  internal  energy  in  the  modified  version 
are  almost  always  below  melt,  the  inclusion  or  exclusion  of  the  SPHASE 
contribution  is  now  primarily  controlled  by  the  tensile  failure  criter- 
ion. If  the  material  density  p is  such  that  p/pq  < AMDM  (pQ  = initial 

density),  the  material  has  failed  and  the  SPHASE  contribution  is  ignored. 


1 

< 

•• 

I 


i 


1 


t 

i 


19.  Von  Holle,  W.  G.  and  Trimble,  J.  J.,  "Shaped  Charge  Temperature 
Measurements,"  Proceeding  of  the  Sixth  Symposium  on  Detonation,  San 
Diego,  August  l3?6. 
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Figure  10.  Cell-Centered  Values  of  the  Jet’s  Specific  Internal  Energy  (MJ/kg)  Computed  by 
the  Modified  Formulation 
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Figure  li.  Modified  Version's  Specific  Internal  Energy  Profiles  from 
the  Stagnation  Region  to  the  Jet  Tip  at  Various  Time 
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The  parameter  AMDM  is  a material  strength  input  value.  Hence,  the 
effects  of  SPHASE  on  a conical  shaped  charge  calculation  are  now  depen- 
dent on  an  input  value.  A sensitivity  study  of  the  calculation  on  this 
parameter  AMDM  should  be  made  with  the  modified  formulation.  However, 
such  a comparison  is  not  included  in  this  report. 

Comparisons  of  the  pressures,  velocities,  specific  total  energies 
and  densities  between  the  two  versions  along  the  axis  of  symmetry  at 
16^s  are  shown  in  Figs.  12-17.  Although  significant  differences  in  the 
above  quantities  are  expected,  especially  within  the  jet,  as  the  result 
of  phase  differences  predicted  by  the  codes  (see  Fig.  9),  the  discrep- 
ancies between  the  versions  are  not  major,  except  for  the  density.  The 
cause  of  the  small  discrepancies  can  be  traced  to  the  pressure  calcula- 
tion. The  effect  of  the  specific  internal  energy  on  the  other  computed 
quantities  is  controlled  by  the  equation  of  state.  In  the  original 
version,  a predetermined  upper  bound  on  the  specific  internal  energy  is 
used  in  the  equation  of  state  for  the  liner  material  in  an  ad  hoc  manner. 
For  shaped-charge  calculations  (code  option  Nl.INER  / 0),  the  maximum 
value  of  the  specific  internal  energy  utilized  in  the  equation  of  state 
for  the  primary  liner  material  is  the  incipient  vaporization  value.  In 
our  example,  the  largest  specific  internal  energy  value  used  in  the 
pressure  calculation  of  the  original  version  is  1.38  MJ/kg,  even  though 
all  the  cells  of  the  jet  have  values  above  1.38  MJ/kg.  In  the  modified 
version  of  HELP,  the  ad  hoc  usage  of  this  upper  bound  is  abandoned. 

Hence,  the  pressure  and  other  variables  are  consistent  with  the  actual 
computed  value  of  the  specific  internal  energy  in  the  modified  version. 

The  pressures  are  compared  in  Fig.  12.  The  graph  corresponding  to 
the  values  computed  by  the  modified  formulation  show  a more  uniform 
decay  to  the  initial  pressure.  The  stagnation  point  overpressure  of 
the  modified  version  is  6%  greater  than  that  of  the  original  version. 

The  stagnation  point  occurs  at  the  same  position  in  both  versions.  The 
graphs  of  the  pressure  and  of  the  radial  velocity  (Fig.  IS)  and  the  com- 
pression (Fig.  17)  correlate  well  for  the  jet  and  slug's  interior. 
Regions  of  high  pressure  correspond  to  large  negative  radial  velocities 
and  values  of  compression  greater  than  one.  The  axial  velocities  of 
the  two  formulations  (Fig.  13)  are  identical  near  the  stagnation  region. 
Both  versions  predict  that  the  physical  inverse  velocity  gradient  is 
still  slightly  present  near  the  jet  tip.  The  large  values  of  the  axial 
velocity  cause  the  slightly  longer  jets  of  the  original  formulation 
(approximately  two  computational  cells).  The  effect  of  the  inverse 
velocity  gradient  on  the  kinetic  energy  can  be  seen  in  Fig.  8.  As  in 
Fig.  13,  the  gradient  is  larger  in  the  modified  version  than  in  the 
original.  The  jet  tip  velocity  values  of  the  original  formulation  (Fig. 
14)  are  somewhat  larger  for  times  greater  than  5us  but  the  percentage 
difference  is  small--at  16ps,  the  difference  is  less  than  2%.  The  rela- 
tive jet  tip  (jet  tip  minus  slug  tail)  velocity  predicted  by  the  modi- 
fied version  at  16|is  is  6.023  km/s.  This  value  is  7.9*0  lower  than  the 
experimental  value  cited  in  Ref.  16.  The  radial  velocities  are  compared 
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Comparison  of  the  Pressures  Along  the  Axis  of  Symmetry  at  16us 


AXIAL  VELOC 


1-igure  13.  Comparison  of  the  Axial  Velocities  Along  the  Axis 
of  Symmetry  at  16ys 


Figure  15.  Comparison  of  the  Radial  Velocities  Along  the  Ajcis  of  Symmetry  at  16us 
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Figure  16.  Comparison  of  the  Specific  Total  Energies  Along  the  Axis 
of  Symmetry  at  16ys 
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in  Fig.  15.  Since  the  radial  velocity  is  modified  so  that  the  jet  ma- 
terial at  zero  overpressure  will  not  expand  radially  (zero  radial  velo- 
cities are  imposed),  an  accurate  assessment  of  this  quantity  is  diffi- 
cult away  from  the  stagnation  region.  However,  in  the  modified  formu- 
lation the  negative  radial  velocities,  which  cause  a compression  of  the 
copper  at  the  axis  of  symmetry,  correlate  well  with  the  nonzero  pressure 
values  in  Fig.  12.  In  the  original  formulation,  this  type  of  correspon- 
dence is  not  as  well  correlated.  The  specific  total  energies  of  the 
two  versions  are  compared  in  Fig.  16.  The  effect  of  the  small  inverse 
velocity  gradient  seen  in  Figs.  8 and  13  is  again  reflected  in  the 
graph  of  the  specific  total  energies.  Fig.  17  shows  a comparison  of 
the  compression  (p/pq)  values.  Throughout  most  of  the  slug  and  jet,  the 

density  values  of  the  modified  code  are  approximately  10%  greater  than 
the  original.  The  density  discontinuity  which  appears  near  the  jet  tip 
in  both  formulations  could  be  numerically  caused.  See  Section  VII.  The 
magnitude  of  the  jet  tip’s  discontinuity  is  larger  in  the  original  ver- 
sion and  so  is  its  associated  undershot.  Before  decreasing  to  the  final 
undershot  value  in  the  modified  version,  the  jet's  density  is  approxi- 
mately 95%  of  the  initial  density  of  the  liner.  No  density  measurements 
have  been  experimentally  determined  for  the  43mm  copper  shaped  charge. 

However,  experimental  results  for  a 36mm  copper  shaped  charge"  pre- 
dict the  density  variation  of  the  unbroken  jet  to  be  from  90%  to  81%  of 
the  initial  liner  density.  Although  the  original  version  seems  to  be  in 
better  agreement  with  Jamet's  results  than  the  modified  version,  a HELP 
simulation  of  the  actual  experiment  with  identical  explosives  and  con- 
figurations would  have  to  be  completed  in  order  to  accurately  ascertain 
which  version  possessed  the  hotter  correspondence.  In  the  stagnation 
region,  the  modified  formulation  gives  improved  valucs-compress ion  (p  > 
p ) of  the  copper.  The  sudden  density  decreases  in  the  slug  region 

predicted  by  the  modified  version  can  be  correlated  to  the  sudden  pres- 
sure decreases.  Sec  Fig.  12. 

Negative  specific  internal  energy  values  are  computed  by  both  for- 
mulations of  the  HELP  code  in  the  uncollapsed  liner  and  along  the  inter- 
face between  the  slug  and  the  explosive.  However,  in  the  modified  for- 
mulation the  number  and  magnitude  of  the  negative  internal  energy  values 
increase.  This  is  consistent,  since  the  original  formulation  added  the 
positive  values  of  the  truncation  error  terms  in  the  specific  kinetic  en- 
ergy calculation  to  the  specific  internal  energy.  Sec  the  discussion  of 
Fig.  8.  The  computed  specific  internal  energy  value  is  the  value  over 
the  initial  specific  internal  energy.  Consequently,  the  actual  internal 
energy  of  a material  in  a computational  cell  may  still  be  positive  even 
though  the  computed  value  is  negative.  Nonetheless,  in  a total  energy 
formulation,  the  computed  negative  specific  internal  energy  values  arc 


20.  Jainet,  F.,  "Measure  de  la  densite  d'un  jet  de  charge  crcuse  en 
curvre  par  radiographic-eclair,"  Saint-l.ouis,  Rapport-Rcrieht  R 101/7<>, 
9. 1. 1970. 
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a consequence  of  the  computed  kinetic  energy  being  greater  than  the  com- 
puted total  energy.  The  numerical  source  of  this  situation  in  the  HELP 
code  must  be  identified  and  corrected  before  any  simulation  can  be  com- 
pletely satisfactory. 


VII.  TEMPERATURE  CALCULATIONS  FOR  THE  43MM  CONICAL  SHAPED  CHARGE 

It  is  not  the  purpose  of  this  section  to  give  a definitive  statement 
on  the  temperature  distribution  within  the  jet  for  a 43mm  shaped  change 
but  rather  to  identify  the  current  problems  in  an  accurate  temperature 
calculation  and  to  present  a first  attempt  in  computing  a temperature 
profile.  The  computed  profiles  do  agree  qual itatively  with  experimental 

■ , 19 

evidence 


Since  the  Tillotson  equation  of  state  used  in  HELP  does  not  compute 

■>] 

temperature  values,  the  temperature  routine  from  the  HULL  code"  was 
utilized.  The  HULL'S  temperature  routine  is  a temperature  fit  which 
uses  data  in  Ref.  22  and  requires  values  of  the  density  and  specific 
internal  energy.  The  coding  used  in  the  HELP  code  is  listed  in  Appen- 
dix C.  This  routine  was  chosen,  since  it  is  used  for  temperature  meas- 
ures in  shaped  charge  simulations  by  the  HULL  code  and  it  is  accessible. 
However,  the  question  of  its  applicability  to  simulations  involving 
phase  transitions  (from  solid,  through  melt,  to  liquid)  is  unanswered. 

A Uefinite  need  does  exist  for  a highly  accurate  copper  equation  of 
state  (see  Ref.  23).  More  sophisticated  equation  of  states,  such  as 
24 

the  BRLGRAY“  equation  of  state,  could  be  used  when  they  arc  incor- 
porated in  the  HELP  code. 

Since  the  temperature  depends  explicitly  on  the  density  and  specific 
internal  energy,  accurate  values  of  these  quantities  arc  essential  to  a 
correct  temperature  calculation.  For  example,  if  the  density  value  is 


21.  klem,  G.  .).,  US  Army  Ballistic  Research  Laboratory,  Private 
Communication,  April  1978. 

22.  Thompson,  S.  L.  and  Lauson,  H.  S. , "Improvements  in  the  Chart  D 
Radiat ion-llydrodvnamic  Code, 1 1 1 : Revised  Analytic  Equation  of  State," 
Sandia  Laboratories  Report  No.  SC-RR-71-0714,  March  1972. 

23.  Edwards,  I..  L.  and  Godfrey,  C.  S. , "Computer  Code  Design  of  Shaped 
Charges,"  Lawrence  Livermore  Laboratory  Report  No.  UCRL-52S98,  October 
1978. 


21.  Wray,  W.  0.  and  Cecil,  R.  A.,  "Modified  Gray:  An  Improved  Three- 
Phase  Equation  of  State  for  Metals,"  HS  Army  Ballistic  Research  Labora- 
tory Contract  Report  No.  299,  April  !97<>.  (AP  #A0252G(1) 
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9.8  x 10  kg/m'  and  specific  internal  energy  value  is  0.46667 
then  a 6.43S  decrease  in  density  and  a 3.  45?o  increase  in  the  specific 
internal  energy  produces  over  a 1S?0  change  in  the  calculated  HULL  tem- 
perature value.  Consider  Tabic  1 which  lists  the  jet's  density  for  the 
modified  version  at  16ps  within  the  computational  cells  from  the  jet  tip 
to  the  stagnation  region.  The  prefix  M with  a number  indicates  that  it. 
is  a free  surface  cell  ( a mixed  cell  composed  of  void  and  copper).  The 
striking  feature  of  the  array  of  numbers  is  that  in  a given  column  the 
density  of  every  free  surface  cell  is  identical.  This  peculiarity 
results  from  the  numerical  algorithm  in  HELP  (see  subroutine  NEWRHO) 
and  is  not  necessarily  a physical  result.  The  same  phenomenon  occurs 
in  bch  the  original  and  modified  formulations.  Thus  computed  tempera- 
ture value  involving  a free  surface  cell  would  be  at  best  suspect. 

Con  -.equently , a comparison  with  Von  liolle  and  Trimbles'  experimental 

25 

je  skin  temperature  “ will  not  be  made.  First,  HELP'S  free  surface 
d'  is  it  y calculation  must  be  analyzed  and  modified,  if  necessary. 

With  an  awareness  of  the  limitations  discussed  above,  we  now  con- 
sider the  temperature  profile  within  the  jet  (pure  cells  only)  with  the 
HULL  temperature  routine.  We  note  that  the  specific  internal  energy 
computed  by  the  HELP  code  must  be  augmented  by  an  initial  value  of 
0.113335748  MG/kg  (corresponds  to  a temperature  of  287. 87K  at  an  initial 
pressure  of  0.101325  MPa)  to  obtain  the  actual  value  of  the  copper's 
specific  internal  energy.  Between  16  and  18.5ns,  the  jet  is  defined 
by  at  most  three  pure  cells  in  the  radial  direction  for  a given  axial 
position.  Figs.  18  and  19  show  the  temperature  profile  at  16  and  18.5 
ys,  respectively,  of  pure  cells  within  the  jet  at  various  radial  posi- 
tions. The  radial  distances  0.26  mm,  0.78  mm  and  1.3  nun  measure  the 
length  from  the  axis  of  symmetry  to  the  center  of  the  first,  second  and 
third  radial  cells,  respectively. 
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The  melting  point  of  copper  under  the  initial  conditions  is 

'j  / 

1356. 6K.“  The  code  predicts  non-ambient  pressures  in  the  jet  along 
most  of  the  axis  of  symmetry  and  ambient  pressures  in  most  of  the 
second  and  third  radial  cells.  Consequently,  the  r = 0.78  mm  and  r = 

1.3  mm  curves  in  both  figures  predict  that  the  jet  is  in  the  solid  state. 
Along  the  axis  of  symmetry  one  cannot  determine  the  actual  stai:e  of  the 
copper  without  more  precise  thermodynamic  data.  Throughout  most  of  the 
l ct , the  results  of  the  modified  HELP  substantiates  Von  Hollo  and 

1 9 

Trimble's  conclusion  of  a solid  jet.  Thus  qualitative  thermal  agree- 
ment between  experimental  evidence  and  computer  simulation  is  achieved 


25~  Von  Holle,  W.  G.  and  Trimble,  .J . .J.,  "Temperature  Measurement  of 
Copper  and  Eutectic  Metal  Shaped  Charge  .lets,"  US  Army  Ballistic 
Research  Laboratory  Report  No.  2004,  1977.  (A0  #B021338L) 

26.  Stull,  D.  R.  and  Prophet,  II.,  (project  directors),  JANAF  Thcrmo- 
chcmical  Tables,  2d  ed,  June  1971. 
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Table  1 . 


Density  Variation  Within  the  .Jet  at  Ions 


Distance  from  Axis  of  Symmetry 


. 26mm 

. 78mm 

1 . 30mm 

1 . 8 2 mm 

M 

12. 132 

M 

12.1 32 

M 

8.9000 

(>  .7517 

M 

8.9000 

6.2893 

M 

8.9000 

6.8630 

M 

8.9000 

7. 1872 

M 

8.9O>'0 

'.5139 

M 

8.9000 

7. "67  2 

M 

8.9000 

7 . 9480 

M 

8.9000 

8.0471 

M 

8.9000 

M 

1 1.312 

8.0855 

8.4719 

M 

11.312 

1 

i 

1 

8.1110 

8.4878 

M 

11.312 

i 

8. 1416 

8.4871 

M 

11.512 

i 

i 

8. j 74  8 

8. 4 "2  5 

M 

11.312 

8.2110 

8.  14  31 

M 

11.512 

i 

8.2516 

8.1116 

M 

1 1.512 

8 . 2899 

8 . 5869 

M 

1 1.312 

i 

8.3201 

8.5754 

M 

1 1.512 

8.3479 

8 . 5634 

M 

11.312 

8.3787 

8.5398 

M 

11.512 

• 

8 . 4203 

8.2"  15 

M 

1 1.312 

M 

9 .511 1> 

8.4490 

8. 1465 

5. 6167 

M 

9 . 5 1 1 (> 

i 

8.4415 

'.9.541 

6.3971 

M 

9. 5116 

j 

8 . 4466 

7.7.157 

7.2789 

M 

9 . 5 1 1 (> 

i 

i 

8.4523 

7.5188 

7. 8684 

M 

9 .511  (> 

i 

i 

8.4590 

7.2S79 

8.2828 

M 

9 .511  (> 

8.4672 

7.  0(i.87 

8. 52', S 

M 

9.51 16 

, 

8.4826 

(i . 8eS7 

8.  (>512 

M 

9.51  10 

j 

8 . 4955 

6.7251 

8 . (>5  55 

M 

9.  51  l<> 
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8 . 5055 

6.6794 

8.  (>595 

M 

9 .511  (> 
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. 26mm 

Di stance 
. 78mm 

from  Axis  of 

1 . 30mm 

.Symmetry 

1 . 82mm 

2 . 34mm 

8.5098 

6.7489 

8 . b662 

M 9.5116 

8.5169 

6.9168 

8.6703 

M 9.5116 

8.5242 

7.1430 

8.6714 

M 9.5116 

8.5315 

7.3787 

8.6683 

M 9.5116 

8.5384 

7.5989 

8.6639 

M 9.5116 

8.5447 

7.7892 

8.6596 

M 9.5116 

8.5503 

7.9460 

8.6555 

M 9.5116 

8.5551 

8.0763 

8.6517 

M 9.5116 

8.5592 

8 . 1 903 

8.6476 

M 9.5116 

8.5633 

8.2976 

8.6428 

M 9.5116 

8.5682 

8.4025 

8.6379 

M 9.5116 

8.5753 

8.4998 

8.6347 

M 9.5116 

8.5862 

8.57.37 

8.6.345 

M 9.5116 

8.6028 

8.6112 

8.6.384 

M 9.5116 

8.6258 

8.6252 

8. 6484 

M 9.5116 

8.6560 

8.6424 

8.6643 

M 9.5116 

8.6986 

8.6752 

8.6847 

M 9.5116 

M 

1 0 . 300 

8.7595 

8. •’211 

8.7095 

8.  1814 

M 

10.300 

8.8457 

8.7933 

8.7785 

8.7366 

M 

10.300 

8.9525 

8.8877 

8.8633 

8.7652 

M 

10.300 

9 . 073.3 

8.9954 

8.9840 

8.8093 

M 

10.300 

9. 1999 

9. 1256 

9. 1073 

8.9486 

9.0024 

9.3208 

9.2644 

9.2399 

9. 1501 

9.1263 

9.4284 

9.5906 

9.3741 

9.3290 

9.2913 

9.5076 

9.4807 

9.4753 

9.4459 

9.4682 

9 . 5560 

9.5382 

9.5370 

9.5223 

9.5073 

9 . 5740 

9.5614 

9.5605 

9. 5528 

9.5553 

I fc  M r t K A 


for  the  first  time.  Finally,  the  temperature  profiles  in  the  jet  de- 
crease with  time  since  the  jet  is  expanding. 
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For  comparison,  the  values  of  the  density  and  specific  internal 
energy  from  the  original  versions  were  used  with  the  HULL  temperature 
routine.  The  corresponding  temperature  values  in  almost  every  pure 
cell  within  the  jet  were  greater  than  3000K.  At  ambient  conditions, 

°6 

the  vaporization  temperature  of  copper  is  2848K.  These  temperature 
results  are  larger  than  those  implied  by  experimental  evidence  by  at 
least  a factor  of  two. 


VIII.  SUMMARY 

We  have  shown  by  a truncation  error  analysis  that  terms  of  the  order 
of  the  truncation  error  in  the  HELP  algorithm  are  included  in  the  kine- 
tic energy  calculation.  In  the  original  HELP  code,  the  updated  values 
of  kinetic  energy  were  computed  as  consequences  of  the  updated  mass  and 
momentum  values.  This  value  is  shown  to  deviate  from  that  computed 
directly  by  a first  order  approximation  of  the  kinetic  energy  equation 
by  first  order  terms  which  depend  quadratical ly  on  the  spatial  deriva- 
tives of  the  velocity,  pressure,  and  elements  of  the  deviator  stress 
tensor.  These  truncation  error  terms  can  severely  alter  the  related 
internal  energy  calculation.  The  effect  of  these  truncation  error  terms 
from  the  kinetic  energy  calculation  on  the  accuracy  of  the  internal  energy 
is  problem-dependent  and  the  criteria  for  such  a determination  is  given  in 
terms  of  the  explicitly  calculated  truncation  error  terms  for  the  specific 
internal  energy.  A method  was  suggested  to  avoid  those  terms  within 
the  confines  of  the  basic  HELP  algorithm.  Although  in  certain  computa- 
tions these  truncation  error  terms  may  remain  negligible,  in  others  they 
can  be  significant  and  produce  spurious  results.  A case  in  point  has  been 
cited  and  illustrated  by  applications  to  problems  in  warhead  mechanics. 

A mesh  refinement  study  for  a copper  wedge  impacting  a perfectly  reflec- 
tive wall  was  made.  The  results  show  tremendous  deviations  as  the  mesh 
size  decreases  in  the  internal  energy  values,  while  other  quantities 
show  relatively  little  deviations.  In  the  modified  formulation,  signi- 
ficantly smaller  deviations  in  the  internal  energy  values  and  conse- 
quently, a more  consistent  deviation  with  the  other  quantities  are 
observed.  Drastic  improvements  in  the  internal  energy  values  are  also 
seen  in  conical  shaped  charge  calculations.  Temperature  profiles  of 
the  jet's  interior  were  computed  using  the  temperature  fit  routine  from 
the  HULL  code.  When  ambient  pressures  are  computed  by  the  code,  the 
corresponding  values  of  the  temperature  are  between  700(K)  and  1100(K) 
which  arc  well  below  that  for  melted  copper.  Thus,  qualitative  thermal 
agreement  is  achieved  for  the  first  time  between  a HF.LP  code  simulation 
and  experimenta1  evidence  for  a conical  shaped  charge. 

The  upper  bound  on  the  specific  internal  energy  (the  value  corres- 
ponding to  the  incipient  vaporization  of  the  material)  in  the  lillotson 
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equation  of  state  for  the  liner  material  has  been  removed  in  the  modi- 
fied formulation.  Now  the  specific  internal  energy  value  used  in  the 
pressure  calculation  is  always  the  value  calculated  by  the  algoritmu. 
Since  the  original  formulation  used  this  ad  hoc  upper  bound,  the  devia- 
tion of  the  pressures,  velocities  and  specific  total  energies  between 
the  two  versions  is  generally  small.  However,  the  density  values  in 
the  modified  formulation  within  the  jet  are  approximately  95%  of  the 
initial  density,  instead  of  approximately  85%  in  the  original  version. 

Although  a comparison  of  a conical  shaped  charge  simulation  with 
and  without  the  SPHASE  portion  of  the  calculation  has  not  been  completed 
with  the  modified  formulation,  the  effects  of  strength  should  be  appre- 
ciable since  the  SPHASE  calculation  is  now  included  by  a significantly 
larger  number  of  pure  cells  (mixed  cells  have  no  strength).  Consequent- 
ly, simulations  performed  for  detailed  jet  information  should  use  the 
strength  option. 

Besides  the  obvious  effects  of  an  improved  internal  energy  calcula- 
tion, another  effect  is  the  increased dependence  on  the  material  con- 
stants used  within  the  code.  A case  in  point  is  in  the  shift  from  the 
internal  energy  dependent  cut-off  value  in  the  SI’HASE  calculation  in 
the  original  formulation  to  a tensile  failure  criterion  in  the  modified 
formulation.  Thus,  the  necessity  of  accurate  material  constants  and 
models  becomes  more  acute. 

As  we  noted  in  Sections  VI  and  VII,  problems  still  exist  within  the 
HELP  code  with  respect  to  an  accurate  internal  energy  and  temperature 
calculation.  The  numerical  problem  related  to  rhe  computation  of  nega- 
tive internal  energies  along  the  slipline  interface  between  the  copper 
and  explosive  as  well  as  in  the  uncol lapsed  liner  must  be  corrected.  The 
determination  of  free  surface  densities  must  also  be  improved.  A veri- 
fied temperature  routine  incorporated  into  the  mollified  version  of  HELP 
is  also  needed.  Once  these  are  accomplished,  the  HELP  code  can  be  used 
with  confidence  to  predict  not  only  the  jet’s  shape  and  iet  tip  velo- 
city, but  also  thermodynamic  quantities  of  the  jet,  such  as  its  density, 
pressure  and  temper. iture. 

In  the  Tl’HASE  section  of  the  original  HELP  algorithm,  the  truncation 
error  terms  were  identified  with  an  explicit  artifical  viscosity  in  the 
internal  energy  calculation.  Consequently,  the  modified  formulation  may 
require  implementation  of  the  artificial  viscosity  option  available  in 
HELP  in  cases  where  the  original  formulation  did  not. 

The  1971  version  of  HELP  lias  a different  energy  formulation.  In 
SPHASE  and  III'IIASE  the  specific  internal  energy  was  updated  directly  by 
a finite  difference  approximation  and  the  specific  total  energy  was 
obtained  as  the  sum  of  the  specific  internal  and  kinetic  energies.  In 
ITHASH,  the  specific  total  energy  was  updated  directly  and  the  specific 
internal  was  obtained  as  the  difference  of  the  specific  totai  and  kine- 
tic energies.  Although  the  proposed  changes  in  this  report  are  no  longer 
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directly  applicable  to  such  a code,  the  basic  concept  developed  in  this 
report  is:  Truncation  error  terms  generated  in  the  calculation  of  a 
kinetic  energy  value  from  the  current  values  of  mass  and  momentum  must 
be  excluded  for  an  acute  overall  energy  calculation. 

A noteworth  feature  of  the  analysis  and  modification  is  that  it  is 
directly  relevant  to  codes  other  than  HELP.  In  fact,  any  code  with  a 
HELP  type  algorithm  can  be  erroneously  affected  by  the  truncation  error 
terms  in  the  kinetic  energy  calculation.  In  particular,  the  same  unphy 
sical  internal  energies  in  the  jet  are  computed  by  the  HULL  code. 

A possible  alternative  solution  to  the  problem  discussed  in  this 
report  is  the  upgrading  of  the  HELP  algorithm  to  second  order  accuracy. 
Such  an  approach  has  proven  beneficial  for  the  standard  particle- in- 
27 

cell  codes.  This  type  of  change  would  require  considerable  analysis 
and  rewriting  of  the  HELP  algorithm.  Present  problems,  which  stem 
primarily  from  the  mixed  cell  treatment  and  which  may  not  be  solved  by 
;•  second  order  algorithm,  should  be  addressed  first. 
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APPENDIX  A 


LISTING  OF  MODIFICATIONS  TO  THE  KINETIC  ENERGY  CALCULATION  IN  HELP 


FfBCEDINa  PAGK  BLAIK-NOT  f ILMKD 


*IDENT  KINENGCOR 
*1  HELPOOM.108 

C ***SDLKEL , SDLKER , SDLKEB , SDLKET  ARE  KINETIC  ENERGIES  OF  A GIVEN 

C MATERIAL  TRANSPORTED  ACROSS  LEFT, RIGHT, BOTTOM  AND  TOP 

C BOUNDARIES  OFAN  INTERFACE  CELL 

C ***DLKEL  IS  KINETIC  ENERGY  OF  MASS  TRANSPORTED  ACROSS  LEFT 

C SIDE  OF  CELL 

*D  HELPQOM. 112 

4 MFGREZ ,REZAMX,REZAIX,REZXMS,REZSIE,REZRHO,SDLKEL, SDLKER, SDLKEB, 

5 f ' KET , DLKEL 
*D  MEMF  ND.45 

3 SFLEFT(4,200) ,SDLKEL (4,200) , SDLKER (4) ,SDLKEB(4) , 

4 SDLKET(4) ,DLKEL(200) 

*1  MEMEXPND.59 

C ***TKBG  IS  SPECIFIC  KINETIC  ENERGY  IN  A CELL 

C ***TKBGM  IS  SPECIFIC  KINETIC  ENERGY  OF  MATERIALS  IN  INTERFACE 

C CELLS 

LEVEL  2,TKEG,TKEGM 

COMMCN/KBGCHK/TKBG (12001)  ,TKEGM(4,800) 

*D  EDIT. 17 

C ***USE  CORRECT  KINETIC  ENERGY 

ESUM=ESUM+XMASS (N ,M) * (SIE (N ,M) +TKE)®1  (N ,M) ) 

*D  EDIT. 20 

C ***USE  CORRECT  KINETIC  ENERGY 

20  ESUM=ESUM+AMX (K) * {TKEG (K) +AIX  (K) ) 

*D  EDIT. 136 
*D  EDIT. 137 

C ***WRITE  CORRECT  KINETIC  ENERGY 

WRITE (KUNITW)  (U (I) ,V(I) ,AMX(I) ,AIX(I) ,TKBG (I) ,P(I) ,MFIAG (I) , 1*1, 

1 KMAX) , (DETIM(I) ,1=1, KDT) 

*D  EDIT. 146 
*D  EDIT. 147 

C ***WRITE  CORRECT  KINETIC  ENERGY 

WRITE  (KUNIIW)  ( (XMASS  (M,L)  ,RHO  (M,  L)  ,SIE  (M,  L)  ,TKEGM  (M,  L) , US  (M,  L)  , 

1 VS(M,L)  ,SAMPY  (M,L)  ,SGAMC  (M,L) , SAMMY  (M,L)  ,SAMMP(M,L)  ,M=1,NMAT)  , 

2 L=1 ,NMXCL£) 

*D  EDIT. 177 

C ***USE  CORRECT  KINETIC  ENERGY 

WSK=WS*TKEG(K) 

*D  EDIT. 190 

C ***USE  CORRECT  KINETIC  ENERGY 

WSX=WS*TKEGM(M,MC) 

*D  EQST.  124 ,EQST.  130 

C DELETE  AD  HOC  UPPER  BOUND  ON  SPECIFIC  INTERNAL  ENERGY 

C IN  THE  PRESSURE  CALCULATION  FOR  THE  LINER  MATERIAL 

*D  ENENCHCK.8 

C ***USE  CORRECT  KINETIC  ENERGY 

EKC=AMX (KK2) *TKEG ( KK2 ) 

*D  ENENCHCK.I3 

C ***USE  CORRECT  KINETIC  ENERGY 


PRICZDINO  BUMUMOT  FILMED 


EKM-XMASS (KLM , MEAN ) ‘TKEGM (KLM, MBAN) 

*1  FILGRD.28 

C “‘OOMPUTE  INITIAL  VALUES  OF  ARRAY  TKEGM 

TKEGM (N,M) »0 . 5* (US (N,M) “2+VS (N ,M) “2) 

*t  FILGRD.32 

C “‘USE  TKFG  AS  A STORAGE  VARIABLE 

TKEG (K) =TKEG (K) +XMASS  (N,M) *TKBGM (N ,M) 

*1  FILGRD.40 

C ‘“COMPUTE  INITIAL  VALUES  OF  ARRAY  TKEG  FOR  MIXED  CELL 

TKEG (K) =TKEG (K) /AMX (K) 

*1  FILGRD.50 

C ‘“SET  MATERIAL  KINETIC  ENERGY  TO  CELL  KINETIC  ENERGY  IN 

C NON-SLIP  CELLS 

TKEGM (N,M)=TKEG(K) 

*1  FILGRD.101 

C “‘COMPUTE  INITIAL  VALUES  OF  ARRAY  TKEG  PGR  PURE  CELLS 

TKEG  (K)  =0 . 5*  (U  (K)  **2+V (K)  “2) 

*D  FILGRD.123 

C “‘USE  CORRECT  KINETIC  ENERGY  TO  COMPUTE  ETH 

162  ETH=ETO+XMASS (N,M) * (SIE (N,M) +TKEGM (N,M) ) 

*D  FILGRD.125 

C ‘“USE  CORRECT  KINETIC  ENERGY  TO  COMPUTE  ETH 

163  ETH=ETH+AMX ( K ) * (AIX (K)  +TKEG (K) ) 

*D  FIGS ET. 51 

C ‘“USE  CORRECT  KINETIC  ENERGY 

ETH=ETH-XMASS (N , M) * ( TKEGM (N , M) +SIE (N , M ) ) 

*D  FIGS ET. 55 

C ‘“USE  CORRECT  KINETIC  ENERGY 

ETH=ETH+XMASS (N , M)  * (TKEGM (N,M) +SIE (N,M) ) 

*D  KPHASE . 168 
*D  HPHASE.169 

C ‘“COMPUTE  CORRECT  CHANGE  IN  KINETIC  ENERGY  BY  MODELLING  THE 

C PARTIAL  DIFFERENTIAL  EQUATION  FOR  THE  CHANGE  OF  KINETIC 

C ENERGY  IN  HPHASE 

DKE=DU*U(K)+DV*V(K) 

C ‘“COMRJTE  CHANGE  IN  SIE  WITH  CORRECT  DKE 

DIE=DE-DKE 
*D  HPHASE. 181 

C ** ‘DELETE  ALTERATIONS  TO  KINETIC  ENERGY  AND  THUS  TO  INTERNAL 

C ENERGY  DUE  TO  CHANGES  IN  THE  VELOCITY  COMPONENTS  IN 

C SUBROUTINE  UVCALC 

257  CONTINUE 

*D  HPHASE. 182, HPHASE. 184 
*D  HPHASE. 191, HPHASE. 196 
*1  HPHASE. 198 

C “‘UPDATE  THE  KINETIC  ENERGY  OF  EACH  MATERIAL  IN  MIXED  CELL 

IF  (XMASS (N,M) .GT.O. ) TKEGM (N,M) =TKEGM (N ,M) +DKE 
‘I  HPHASE. 200 

C ‘“UPDATE  KINETIC  ENERGY  ARRAY 

TKEG (K) =TKEG (K) +DKE 

<>.S 


i 

I 


*D  INPUT. 163 
*D  INPUT. 164 

C ***READ  THE  CORRECT  KINETIC  ENERGY 

READ  (KUNITR)  (U(I) ,V(I) ,AMX(I) ,AIX(I) ,TKEG(I) ,P(I) ,MFLAG(I) ,1=1, 

1 KMAX) , (DETIM(I) ,I=1,KDT) 

*D  INPUT. 173 
*D  INPUT’.  174 

C ***READ  CORRECT  KINETIC  ENERGY 

READ  (KUNITR)  ( (XMASS  (M,L)  ,RHO(M,L)  ,SIE  (M,L)  ,TKEGM(M,L) , US (M,L) , 

1 VS(M,L)  ,SAMPY(M,L)  ,SGAMC  (M,L) , SAMMY (M,L)  ,SAMMP{M,L)  ,M«1,IMAT) , 

2 L=1,M4XCLS) 

*1  NEWFIG.16 

C ***SET  KINETIC  ENERGY  OF  CELL  NOT  IN  USE  TO  ZERO 

TKEGM (NN ,M) “0 . 0 
*1  NEWMIX.14 

C ***DEFINE  KINETIC  ENERGY  FOR  NEW  MIXED  CELL 

TKEGM (MO, M)=TKBG(Kj 
*D  RNDOFF . 27 

C ***USE  THE  CORRECT  KINETIC  ENERGY 

OUTK=DM*TKEGM (N ,M) 

*D  SETUP.  136 
*D  SETUP. 137 

C ***WRITE  THE  CORRECT  KINETIC  ENERGY 

WRITE (KUNITW)  (U(I) ,V(I) ,AMX(I) ,AIX(I) ,TKEG (I) ,P(I) ,MFIAG (I ) , 1=1 , 

1 KMAX) , (DETIM(I) ,I=1,KDT) 

*D  SETUP. 147 
*D  SETUP. 148 

C ***WRITE  THE  CORRECT  KINETIC  ENERGY 

WRITE  (KUNITW)  ' (XMASS  (M,L)  ,RHO(M,L)  ,SIE(M,L)  ,TKEGM(M,L)  ,US(M,L) , 

1 VS  (M,L)  .SAMPY (M,L)  ,SGAMC(M,L) , SAMMY (M,L)  ,SAMMP(M,L)  ,M=1,NMAT)  . 

2 L=1,NMXCLS) 

*1  SETUPA.112 

C *** INITIALIZE  KINETIC  ENERGY  MATERIAL  ARRAY 

TKEGM (N,M) =0.0 
*1  SETUPA.  123 

C ***INITIALIZE  KINETIC  ENERGY  GRID  ARRAY 

TKBG(K)=0.0 
*D  SPHASE . 628 

C ***COMPUTE  CORRECT  CHANGE  IN  KINETIC  ENERGY  BY  MODELLING  THE 

C PARTIAL  DIFFERENTIAL  EQUATION  FOR  THE  CHANGE  OF  KINETIC 

C ENERGY  IN  SPHASE 

DKE=DELU*UKT+DELV*VKT 
*D  SPHASE. 629 

C ***COMKJTE  CHANGE  IN  SIE  WITH  CORRECT  DKE 

DIE=DE-DKE 
*1  SPHASE. 644 

C ***UPDATE  THE  SPECIFIC  KINETIC  ENERGY  OF  EACH  MATERIAL  IN 

C A MIXED  CELL 

TKEGM (MM , MKF) -TKEGM (f« , MKF ) +DKE 
*1  SPHASE. 654 


A.  •* 


C “‘UPDATE  KINETIC  ENERGY  ARRAY 

TKEG (K) “TKEG (K) +DKE 
‘D  SPHASE . 690 
*D  SPHASE. 691 

C ‘“USE  CORRECT  KINETIC  ENERGY 

E=E+XMASS (NN2,MM2) * (TKEGM(NN2,MM2) +SIE (NN2,MM2) ) 

‘D  SPHASE. 694 

C ‘“USE  CORRECT  KINETIC  ENERGY 

560  E=E+AMX (KK2 ) * (TKEG (KK2 ) +AIX (KK2 ) ) 

*1  TPHASE. 23 

C “‘SET  KINETIC  ENERGY  FLUX  OF  EACH  MATERIAL  AT  AXIS  TO  ZERO 

SDLKEL(N,J) =0. 

*1  TPHASE.61 

C ‘“SET  KINETIC  ENERGY  FLUX  OF  CELL  AT  AXIS  TO  ZERO 

DLKEL (J)=0. 

*1  TPHASE . 82 

C ‘‘‘INITIALIZE  KINETIC  ENERGY  FLUX  OF  EACH  MATERIAL  AT  BOTTOM 

C TO  ZERO 

SDLKEB (N)=0. 

*D  TPHASE. 100 

C “‘USE  CORRECT  KINETIC  ENERGY  FOR  BOTTOM  FLUX  OF  PURE  CELL 

DLKEB=TKEG (K) *AMMY 
DELEBK=DLKEB 
*1  TPHASE. 105 

C “‘INITIALIZE  KINETIC  ENERGY  FLUX  AT  BOTTOM 

DLKEB=0 . 

*1  TPHASE. 11 3 

C “‘USE  CORRECT  KINETIC  ENERGY  FOR  BOTTOM  FLUX  OF  EACH  MATERIAL 

SDLKEB  (N)  =TKE)GM (N,M)  ‘SAMMY  (N,M) 

*D  TPHASE. 11 4 

C “‘USE  CORRECT  KINETIC  ENERGY  FLUX  OK  EACH  MATERIAL  TO 

C DETERMINE  FLUX  OF  WHOLE  CELLFOR  KINETIC  AND  TOTAL 

C ENERGY  FLUXES 

DLKEB =DLKEB  +SDLK  EB (N) 

DELEB=DELEB+SDELEB (N) +SDLKEB (N) 

*D  TPHASE. 116 

C “‘USE  CORRECT  KINETIC  ENERGY  FOR  OUTFLCW  OF  KINETIC  ENERGY 

C AT  BOTTOM  OF  GRID 

OUTKE  (N)  =OUTKE  (N)  -SDLKEB  (N) 

*1  TPHASE. 149 

C ‘“SET  CORRECT  KINETIC  ENERGY  FLUX  AT  BOTTOM  TO  ZERO 

DLKEB=0. 

*1  TPHASE.  158 

C ‘“INITIALIZE  SPECIFIC  KINETIC  ENERGY  FLUXES  AT  TOP 

C AND  RIGHT  TO  ZERO 

EKAMPY=0.0 
EKAMMP=0 . 0 
*D  TPHASE. 218 

C ‘“USE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  KINETIC 

C ENERGY  FLUX 


7U 


EKAMPY=TKEGM (MFK ,MFL-10G) 

C ***USE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  TOTAL 

C ENERGY  FLUX 

EAMPY=SIE (MFK,MFL-100) +EKAMPY 
*D  TPHASE. 227 

C ***USE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  KINETIC 

C ENERGY  FLUX 

330  EKAMPY-TKEG (M) 

C ‘“USE  CORRECT  KINETIC  ENERGY  TO  DETFRMINE  SPECIFIC  TOTAL 

C ENERGY  FLUX 

DAMPY =AIX  (M)  +EKAMPY 
*1  TPHASE. 236 

C ‘“USE  CORRECT  KINETIC  ENERGY  FLUX  TO  FIND  FINAL  VALUE  AT  TOP 

DLKET-AMPY*EKAMPY 
*D  TPHASE. 230 

C ‘“USE  CORRECT  KINETIC  ENERGY  TOR  OUTFLOW  OF  KINETIC  ENERGY 

C AT  TOP  OF  GRID 

OUTKE  (MFK)  TOUTKF  (MFK)  +AMPY*TKEG  (M) 

*D  TPHASE. 293 

C “ *'JSE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  KINETIC 

C ENERGY  FLUX  AT  RIGHT  FOR  MIXED  CELL 

EKAMMPOTKEGM (MFK ,MFI  -100) 

C ‘“USE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  TOTAL 

C ENERGY  FLUX  FOR  MIXED  CELL 

E.AMMP  =S I E (MFK , MFL- 100 ) + EKAMMP 
*D  TPHASE. 30 i 

C ‘“USE  CORRECT  KINETIC  ENERGY'  TO  DETERMINE  SPECIFY  KINETIC 

C ENERGY  FLUX  AT  RIGHT  FOR  PURE  CELL 

460  EKAMMP=TKEG (M) 

C ‘“USE  CORRECT  KINETIC  ENERGY  TO  DETERMINE  SPECIFIC  TOTAL 

C ENERGY  FLUX  FOR  PURE  CELL 

EAMM'3=AI  X (M)  +EKAMMP 
*1  'TPHASE.  311 

C “‘USE  CORRECT  KINETIC  ENERGY  FITJX  TO  FIND  FINAL  VALUE  AT  RT. 

DI.KER=AMMP*EKAMMP 
*D  TPHASE. 314 

C “‘USE  C <RECT  KINETIC  ENERGY  FOR  OUTFLOW  OF  KINETIC  ENERGY  AT 

C R1GI"  IF  GRID 

OUTKE (MFK) =OUTKE (MFK) +AMMP*TKEG (M) 

*1  TPHASE. 319 

C ‘‘‘INITIALIZE  TOP  AND  RIGHT  FLUXES  OF  KINETIC  ENERGY  TO  ZERO' 

DLKKT-n). 

DLKER=0 . 

*1  TPHASE. 341 

C ** ‘USE  CORRECT  KINETIC  ENERGY  FOR  KINETIC  ENERGY  FLUX 

EKAM I !Y  - IK  IX  , ( KT ) 

*1  TPHASE. 344 

C * * ‘USE  CORRECT'  KINETIC  ENERGY  FOR  KINETIC  ENERGY  FIUX 

EKAMPY=TKEGM  (N  ,MT) 

*!  TPHASE . 356 


c 


“‘USE  CORRECT  KINETIC  ENERGY  FOR  KINETIC  ENERGY  FLUX 
EKAMMP=TKEG (KR) 

*1  TPHASE . 359 

C ***USE  CORRECT  KINETIC  ENERGY  FOR  KINETIC  ENERGY’  FLUX 

EKAMMP=TKEGM (N ,MR) 

*D  TPHASE. 362 
*D  TPHASE. 363 

***USE  CORRECT  SPECIFIC  KINETIC  ENERGY  FLUXES 
560  WSA=EKAMFY 
WSB=EKAMMP 
*1  TPHASE. 365 

C ‘“INITIALIZE  MIXED  CELL  TRANSPORT  VARIABLES 

SDLKET (N) =0 . 

SDLKER (N) =U . 

*D  TPHASF . 384 

C ‘“COMPUTE  CORRECT  KINETIC  ENERGY  FLUX  FOR  EACH  MATERIAL 

SDLKET (N) =SAMPY (N,MK) *WSA 
C ‘“COMPUTE  CORRECT  TOT.AL  ENERGY  FLUX 

DELET=DELET+SDELET (N ) +SDLKET  (N) 

C “‘COMPUTE  CORRECT  KINET’IC  ENERGY  FLUX  OF  CELL 

DLKET=DLKET+SDLKET  (N ) 

*D  TPHASE. 397 

C “‘COMPUTE  CORRECT  KINETIC  ENERGY  FLUX  FOR  EACH  MATERIAL 

SDLKER (N ) =SAMM  P (N ,MK) *WSB 
C ‘“COMPUTE  CORRECT  TOTAL  ENERGY  FU.JX 

DELER-DELER+SDELER (N) +SDLKER (N) 

C ‘“COMPUTE  CORRECT  KINETIC  ENERGY  FU;X  OF  CELL 

DLKER=DLKER + SDLKER (N) 

*D  LCHANGE.I6 
*D  TPHASE. 4 30 

C ‘“USE  CORRECT  KINETIC  ENEiAJY  FAIR  OUTFLOW  OF  KINETIC  ENERGY  AT 

C TOP  OF  GRID 

620  OUTKE  (N)  =OUrrKE  (N)  + SDLKET  (N) 

*1  TPHASE. 436 

C “‘SET  KINETIC  ENERGY'  FLUX  ’IO  ZERO 

DLKCT=0 . 0 
*D  LCHANGE . 1 7 
*D  TPHASE  .4  57 

C ‘“USE  CORRECT  KINETIC  ENERGY  FOR  OUTFLOW  OF  KINETIC  ENERGY  AT 

C BOTTOM  OF  GRID 

660  OUTKE (N) =OUTKE (N) +SDLKEB (N) 

*1  TPHASE. 4 62 

C ‘“SET  KINETIC  ENERGY  FUJX  'IT'  ZERO 

DIXER--0.0 
*D  TPHASE. 465 

C ‘“USE  CORRECT  KINETIC  ENERGY 

WSA-TKBG (K) 

*1  TPHASE. 467 

C ‘“COMPUTE  TOTAL  CHANGE  IN  KINETIC  ENERGY  IT.'R  PURE  CELL 

S IGKE=-DLKET-DLKER+DIJ<EB+DLKEL  (J ) 


n n 


*1  TPHASE.47C 

C ***INITIALIZE  VALUE  OF  EKNEW  TO  ZERO 

EKNEW=0 . 

*1  TPHASE .473 

C ‘“COMPUTE  NEW  VALUE  OF  KINETIC  ENERGY  FOR  PURE  CELL 

EKNEW= (SIGKE+AMX (K) *TKEG (K) )/WS 
*1  TPHASE. 486 

C ‘“OOMRJTE  NEW  VALUE  OF  KINETIC  ENERGY  FOR  MIXED  CELL 

SSIGK£=-SDI KET (N) -SDLKER (N) +SDLKEL (N, J) +SDLKEB (N) 

*1  TPHASE .492 

C ‘“SET  KINETIC  ENERGY  OF  CELL  TO  ZERO 

TKEGM (N,MK)=0. 

*1  TPHASE. 500 

C ‘“TRANSFER  CORRECT  KINETIC  ENERGY  VALUE 

EKNW^EKNEW 
*1  TPHASE. 503 

C ‘“COMPUTE  TOTAL  CHANGE  IN  KINETIC  ENERGY  FOR  MIXED  CELL 

EKNW=  (SSIGKE+XMASS  (N ,MK)  ‘TKEGM  (N,MK) ) /SWS 
*D  TPHASE. 504 ,TPHASE. 516 

C ‘“LINES  504-516  WERE  DELETED  SINCE  KINETIC  ENERGY  EQUATION 

IS  BEING  MODELLED  AND  THERMALIZED  ENERGY  IS  A PROPERTY  OF 
ORIGINAL  SCHEME 

720  SWSA=(SIE (N,MK)*XMASS (N,MK) +SWSB)/SWS 
*1  TPHASE. 52 5 

C ‘“SET  COMPUTED  KINETIC  ENERGY  TO  ARRAY  VALUE  FOR  A MIXED  CELL 

TKEGM (N,MK)=EKNW 
*D  TPHASE. 538 

C ‘“COMPUTE  NEW  VALUE  OF  INTERNAL  ENERGY  BASED  ON  CORRECT  K.E. 

S IENEW= (AIX (K) *AMX (K ) +WSB-SIGKE ) /WS 
‘I  TPHASE. 547 

C ‘“SET  COMPUTED  KINETIC  ENERGY  TO  ARRAY  VALUE  FOR  A PURE  CELL 

TKEG (K)=EKNEW 
*D  TPHASE. 569 

C ‘“COMPUTE  TOTAL  ENERGY  IN  GRID  WITH  CORRECT  KINETIC  ENERGY 

ENEPGY=ENERGY+AMX  (I_JD)  * (AIX  (LJD)  +TKEG  (LTD) ) 

*D  TPHASE. 564 
*D  TPHASE. 56 5 

C ‘“COMPUTE  TCTAI.  ENERGY  IN  GRID  WITH  CORRECT  KINETIC  ENERGY 

ENERGY=ENERGY+XMASS(NJ,MJ)*(SIE(NJ,MJ)+TKEC;M(NJ,MJ)) 

*1  TPHASE. 57 3 

C '“WRITE  KINETIC  ENERGY'  FLUXES  FOR  CELL  CENTERED  QUANTITIES 

WRITE  (6,1330)  DLKCT , DLKER , DLKEB , DLKEL ( J ) 

*1  TPHASE. 582 

C ‘“WRITE  KINETIC  ENERGY  FLUXES  FOR  EACH  MATERIAL  IN  MIXED  CELL 

WRITE  (6,) 340)  (N ,SDi.KET (N) , SDLKER (N) ,SDLKEB (N) ,SDLKEL (N, J ) ,N=1 , NM 
LAT) 

*1  1PHASE . 588 

C ‘“SET  KINETIC  ENERGY  FLUXES  TO  ZERO  FOR  EACH  MATERIAL 

SDLKEB (N) =0 . 

SDLKEL (N, J) =0 . 


*1  TPHASE. 608 

C ***STORE  FLUXES  TO  BE  USED  IN  TRANSPORT  OF  KINETIC  ENERGY 

SDLKEB  (N)  =SDIJ<ET  (N) 

SDLKEL (N,  J ) =SDLKER (N) 

*1  TPHASE. 626 

C ***STORE  FLUXES  FOR  KINETIC  AND  INTERNAL  ENERGIES  WHICH  USE 

C CORRECT  KINETIC  ENERGY 

SDLKEB (MFK) = E KAMPY  * AMPY 
SDELEB (MFK) =DELET-SDLKEB (MFK) 

*1  TPHASE. 633 

C ***’ STORE  FLUXES  FOR  KINETIC  AND  INTERNAL  ENERGIES  WHICH  USE 

C CORRECT  KINETIC  ENERGY 

SDLKEL (MFK ,J) =EKAMMP*AMMP 
SS IGC (MFK , J ) =DELER-SDLKEL (MFK , J ) 

*1  TPHASE. 643 

C ***STORE  KINETIC  ENERGY  FLUX 

DLKEL (J) =DLKER 
*1  TPHASE. 647 

C ***STORE  KINETIC  ENERGY  FLUX 

DLKEB=DLKET 
*D  TPHASE. 698 

C ***USE  CORRECT  KINETIC  ENERGY 

WS=TKEG  (K) 

*J  TPHASE. 71 5 

C ***INITIALIZE  KINETIC  ENERGY  SUM  TO  ZERO 

EKSUM=0 . 

*D  TPHASE. 728 

C ***USE  CORRECT  KINETIC  ENERGY 

WS-TKEGM  i'N  , M ) 

*1  TPHASE. 740 

C ***SET  KINETIC  ENERGY  OF  MIXED  CELI.  TO  ZERO 

TKDGM (N ,M) =0 . 

*1  TPHASE . 744 

C ***SUM  KINETIC  ENERGIES 

EKSUM=EKSUM+TKEGM(N,M) *XMASS(N,M) 

*1  TPHASE. 753 

C ***DEF INF,  KINETIC  ENERGY  OF  CELI, 

TKEG  (K)  =EKSUM/XMSUM 
*1  TPHASE. 758 

C ***SET  KINETIC  ENERGY  OF  CELL  TO  ZERO 

TKBG (K) =0 . 

*1  TPHASE. 771 

C ** *SET  KINETIC  ENERGY  OF  MIXED  CEIJ,  TO  ZERO 

TKBGM  (N  ,MK)  =0  . 

*1  TPHASE. 78 i 

C ***SET  KINETIC  HIERGY  OF  CELI,  TO  ZERO 

TKEG  (K)  =0. 

*J  TPHASE . 857 

1330  FORMAT  (7ii  DLKET= IPEI5.8, 6X, 6HDLKER=lPEl 5 .8, 6X , 6HDLKEB=lPEl 5. 8, 10(1 
I DLKEL (J) -1PE15. 8) 


I 


1340  FORMAT  (13, 8H  SDLKFT=1PE15.8,6X,7HSDLKER=1PE15.8,6X,7HSDLKEB=1PE15 
1.8, 6X, 7HSDLKEL=lPE15 .8) 

*D  UVMOD.69 

C ***DELETE  ALTERATIONS  TO  KINETIC  ENERGY  AND  THUS  TO  INTERNAL 

C ENERGY  DUE  TO  CHANGES  IN  THE  VELOCITY  COMPONENTS  IN 

C SUBROUTINE  UVMCO 

C ***LINES  69,70,92,83,85,87-91  OF  UVMOD  WERE  DELETED 

*D  UVMOD. 70 
*D  UVMOD. 82 
*D  UVMOD. 83 
*D  UVMOD. 85 
40  CONTINUE 
*D  UVMOD. 87, UVMOD. 91 


*ID  CORMAP 
*D  MAP. 13 

DATA  XUM  /2H  . , 2H— , 2H-A, 2H-B, 2H-C, 2H-D, 2H-E, 2H-F, 2H-G, 2H-H, 2H-I , 2 


r>  o 


mm 


*ID  REZJAS 
*D  MAIN. 23 

REZONE  AT  CYCLE  REYC9 
USER  MUST  SUPPLY  REYC9  IN  NEXT  STATEMENT 
IF(NC  -NE.  REYC9)  GO  TO  30 
* I DENT  CHREZ 

*D  REZONE. 159, REZONE. 160 

C SET  NJMAX=JEXTY  WHERE  JEXTY  = THE  NUMBER  OF  ROWS  TO  BE  REZONED 
NJMAX=JEXTY 

IF  {MOD (JEXTY, 2)  .EQ.  0)  GO  TO  210 
*D  REZONE. 172 

210  NJMAX1 =NJMAX/ 2 
NJMAX3=NJMAX+1 
*D  REZONE. 182 

IF  (I.GE.  NJMAX)  GO  TO  225 

IF  (I  .LT.  NJMAX  .AND.  IVARDY  .EQ.  1)  GO  TO  220 
*1  REZCNE. 188 

GO  TO  230 

225  TY (N,M) =TY (N,M) -NJMAX1 
*D  REZONE. 196 

IF (I  .GE.  NJMAX)  GO  TO  255 
IF (I . LT.  NJMAX  .AND.  IVARDY. EQ.  1)  GO  TO  250 
*1  REZONE. 202 

255  YP (M) =YP (M) -NJMAXl 
*D  REZONE. 238 

DO  320  J=l, NJMAXl 
*1  REZCNE. 244 

NJMAX2=NJMAX1+1 
DO  324  J =N JMAX3 , JMAX 
DO  325  I=1,IMAX 
K= (NJMAX2-1) *IMAX+I+1 
M=(J-1) *IMAX+I+1 
MFLAG (K) =MFLAG (M) 

AIX(K) =AIX(M) 

U (K)  =U  (M) 

V (K)  =V  (M) 

TKEG  (K)  =TKEG  (M) 

AMX  (K)  =AMX  (M) 

IF (CYCFH3  .LT.  0)  GO  TO  326 
STRSRR(K)  =STRSRR{M) 

STRSRZ (K) -STRSRZ (M) 

STRSZZ (K) =STRSZZ (M) 

326  IF (IDRT  .LE.  0)  GO  TO  327 
DETIM (K) =DETIM (M) 

C **  ZERO  OUT  ORGINAL  CELL  *** 

327  AIX (M) =0 . 0 
AMX (M) =0.0 
U(M) =0.0 

V (M) =0.0 
TKEG (M) =0.0 


PRECEDING  PACE  BUNK-NOT  FIIMED 


IFUDRT  .LE.  0)  00  TO  328 
DETIM(M) =0.0 

328  IF  (CYCFA3  .LT.  0)  GO  TO  329 
STRSRR(M) =0.0 

STRSR2 (M) =0.0 
STRSZZ (M) =0 . 0 

329  CONTINUE 
325  CONTINUE 

NJMAX2=NJMAX2+1 
324  CONTINUE 
*D  REZONE. 248 

DO  330  J=1 , NJMAX1 
*D  REZONE. 250, REZONE -251 
NJMAX2=NJMAX1+1 
NJMAXS=NJMAX 
NJMAX4 =JMAX~N JMAX 1 
DO  340  J=HJMAX2 , NJMAX4 
Y ( J ) =Y  (NJMAXS-i  i ) 

340  NJMAXS=NJMAXS+1 
NJMAX5=NJMAX4+1 
YDELT=Y (NJMAX4 ) -Y (NJMAX4-1) 

DO  345  J=NJMAX5,JMAX 
345  Y (J) =Y (J-l ) +YDELT 
*1  REZONE. 257 

DO  370  J=N JMAX 5 , JMAX 
*D  REZONE. 259 
*D  REZONE. 282 

DO  380  J=NJMAX5 , JMAX 
*D  REZCNE. 28 5, REZONE. 287 
I2=NJMAX5+2 

IF  (JDBT  .GT.  0 ) JDBT=NJMAX5+1 
IF  (JTYF  .GT.  0)  JDTP=NJMAX5+1 
*D  REZONE. 303 
1TY=, 15,/) 


* I DENT  CORREZ 
*D  MEMEXPND. 58 
3 

*1  REZONE. 133 

TKEG (K) =0 . 0 
*1  REZONE. 141 


REZRHO{4, 200) ,TTKM(2,4) 


TKEGM(N,MFK)=0.0 
*7  REZONE. 265 

TKEG (K) =0.0 
*1  REZONE. 274 


TKBGM (N,MFK) =0 . 0 


M, 


*D  OOMPRS . 26 , COMPRS . 30 
*B  OOMPRS.  31 

AIX (K) = (AIX (L) *AMX (L) +AIX (M) *AMX (M) ) /WSA 
TKEG  (K)  = (TKEG  (L)  *AMX  (L)  +TKEG  (M)  *AMX  (M) ) /W SA 
*1  COMPRS. 48 

TKEG  (M)  =0 . 0 
*1  COMPRS.  59 

TKEG (L) =0.0 
*1  COMPRS. 80 

TTKM(II,JJ)=0.0 
*1  COMPRS. 103 

TTKM  (JJ , II ) =TKEGM  (II  ,N) 

*1  OOMPRS. 109 

TKEGM(II,N)=0.0 
*1  COMPRS . 126 

TTKM(1,MFM)=TKEG(M) 

*1  COMPRS. 138 

TTKM (2,MFL)  =TKEG  (L) 

*1  COMPRS. 152 

TKEG (K) =0 . 0 

*D  COMPRS. 16 5, COMPRS. 170 
*B  OOMPRS . 17 1 

SIE (II ,N) = (SSIE (1,11) *XMAS (1 , II) +SSIE (2,11) *XMAS (2,11))  /WSB 
TK EGM ( 1 1 , N ) = ( TTKM ( 1 , 1 1 ) * XMAS ( 1 , 1 1 ) +TTKM ( 2 , 1 1) * XMAS ( 2 , 1 1 ) ) /WSB 
*1  OOMPRS. 173 

TKEG (K) =TKEG  (K)  +TXEGM (II ,N) *XMASS (II ,N) 

*1  OOMPRS. 177 

TKEG (K) =TKHG (K)  /WSA 
*D  OOMPRS. 184 
*D  OOMPRS. 185 
*1  COMPRS.  181 
TKE=0 . 0 
*1  COMPRS. 186 

TKE=TKE+TKEGM(II ,N) *XMASS (II ,N) 

*1  COMPRS. 191 

TKEG (K) =TKE/WSA 
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The  HULL  temperature  routine  uses  the  absolute  specific  internal 
energy  value  and  not  the  specific  internal  energy  value  above  the  ini- 
tial value.  Consequently,  this  routine  augments  the  specific  internal 

energy  value  from  the  HELP  code  by  0.113335748  MJ/kg  which  cor- 
responds to  a temperature  of  287. 87K  at  ambient  pressure  of  0.101325 

3 3 

MPa  and  ambient  density  of  8.9  x 10'  kg/m".  The  routine  has  two  pur- 
poses: to  calculate  a temperature  value  for  the  copper  material  accord- 
ing to  the  HULL  temperature  routine  and  to  list  the  cell  values  of  cer- 
tain quantities  associated  with  the  copper  material.  It  was  appended 
to  the  MAIN  subroutine  by  an  UPDATE  directive  INSERT. 


*ID  KJLLTEMP 
*1  MAIN.  16 

WRITE (6, 1025) 

1025  FORMAT ( 1H1 , 5X , 1HI , 3X , 1HJ , 6X , 4HMASS , 9X , 5RU-VEL , 8X , 5HV-VEL , 9X , 3HSI E , 
* 10X , 3HSKE , 8X , 7HU2+V2/2 , 6X , 7HDENSITY , 7X , 4HTEMP , 1 OX , 3HSTE) 

REFRHO=8 . 9 
Fl=-1285. 

F2=-10000. 

F3=2.54E-07 
F4=1.63E-07 
EMP=1.41E10 
ES=5. 31E10 
DO  1021  1=1 ,IMAX 
DO  1021  J=1 , JMAX 
K=(J-1) *IMAX+I+1 

IF  (MFLAG(K)  .GT.  100)  GO  TO  1023 
IF  (MFLAG(X)  .NE.l)  GO  TO  1021 
C PURE  CELL 

SKE=0 . 5* (U (K) **2+V (K) **2) 

E=AIX (K) 

DENSITY=AMX (K) /(TAU (I) *DY (J) ) 

KEY=1 

GO  TO  1031 
1032  CONTINUE 

TENG=AIX  (K)  +TKEG  (K) 

WRITE (6, 1022) I, J,  AMX(K) ,U(K) ,V(K) ,AIX(K) ,TKEC (K) ,SKE, 

* DENSITY, TEMP, TENG 

1022  FORMAT ( 1H  ,1HP,  2{1X,I4) ,9  (1X,1PE12.4) ) 

GO  TO  1021 

C MIXED  CELL 

1023  M=MFIAG  (K) -100 

IF (SIE ( 1 ,M)  .EQ.  0.0)  GO  TO  1021 
N=1 

SKE=0 . 5* (US (N ,M) **2+VS (N,M) **2) 

DENSITY-RHO(l,M) 
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E=SIE(1,M) 

KEY=2 

GO  TO  1031 
1033  CONTINUE 

1027  TENG=SIE (N,M)+TKEGM (N,M) 

WRITE  (6,1024)  I,J,  XMASS(N,M)  ,US(N,M)  ,VS(N,M)  ,SIE(N,M)  ,TKEGM(N,M 
* ) ,SKE,  DENSITY, TEMP ,TENG 
1024  FORMAT  (1H  ,1HM,  2 (1X,I4) ,9 (IX, 1PE12.4)  ) 

GO  TO  1021 
1031  CONTINUE 

C HULL  TEMPERATURE  ROUTINE 
E=E+1 . 13335748E9 
RMU-DENSITY/REFRHO-1 
RMU=AMAXl (RMU,-1 .0) 

RMU=AMIN1 (RMU, 0.679) 

1 S=AMAX1  ( ( E-EXP)  *1 . E-03 ,0.0) 

j EE=AMIN1 (E,EWP+AMAX1 (E-EMP-ES,S) ) 

1 IF (RMU. LT. 0.0)  GOTO  100 

1 CONl= (F2*RMU+F1 ) *RMU 

CON2=F3 

i GOTO  200 

100  CON1=0. 0 

CON2~F3+F4*RMU 
200  TEMP=CONl+CON2*EE 
DTDE=CON2 

IF(E.CT.EMP.AND.E.LT. (EMP+ES) ) DTDE=DTDE*1 . OE-03 
IF (TEMP  .GE.10.0)  GO  TC  11 
IF(OONl.EQ.O.O)  GOTO  11 
CON2=10 . 0*CON2/ ( 10 . 0-CONI ) 

TEMP=CON2*E 
DTDE=CON2 
11  CONTINUE 

IF  ( KEY  .EQ.  1)  GOTO  1032 
GO  TO  1033 
1021  CONTINUE 


LIST  OF  SYMBOLS 


- specific  kinetic  energy  (J/kg) 

- mass  per  unit  length  (kg/m) 

- time  (s) 

- x-component  of  velocity  (m/s) 

- y-component  of  velocity  (m/s) 

- Cartesian  spatial  coordinate  (m) 

- Cartesian  spatial  coordinate  (m) 

2 

- cross  sectional  area  of  computational  cell  (m  ) 

- boundary  of  computational  cell  traversed  in  a positive  sense  (m) 

- specific  total  energy  (.l./kg) 

- specific  internal  energy  (J/kg) 

- length  of  cell  boundary  (m) 

- order  symbol 

- pressure  (Pa) 

- stress  deviator  tensor  (Pa) 

- transported  mass  per  unit  length  (kg/m) 

- O.SpuAx 

- 0.5pu2At 

3 

- density  (kg/m  ) 

- time  increment  (s) 

- grid  increment  in  the  x-direction  (m) 

- grid  increment  in  the  y-direction  (m) 


SUBSCRIPTS 


denotes  the  number  of  the  cell  in  the  x-Jirection 
denotes  the  number  of  the  cell  in  the  y-direction 
denotes  the  number  of  time  increments 

SUPERSCRIPTS 

denotes  the  computational  cell  boundary  above  its  center 

denotes  the  computational  cell  boundary  below  its  center 

denotes  generic  computational  cell  boundary 

denotes  the  computational  cell  boundary  to  the  left  of  its  center 

denotes  the  computat ional  cell  boundary  to  the  right  of  its 

center 

denotes  value  of  a variable  at  the  completion  of  S PHASE 
denotes  value  of  a variable  at  the  completion  of  Hl'IIASH 
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